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P-] ■ ABSTRACT 

-d ■ 

We present models for the formation of terrestrial planets, and the collisional evolution of debris disks, in planetary systems that con- 
tain multiple marginally unstable gas giants. We previously showed that in such systems, the dynamics of the giant planets introduces 
a correlation between the presence of terrestrial planets and cold dust, i.e., debris disks, which is particularly pronounced at A ~ 70/um. 
Here we present new simulations that show that this connection is qualitatively robust to a range of parameters: the mass distribution 
of the giant planets, the width and mass distribution of the outer planetesimal disk, and the presence of gas in the disk when the giant 
planets become unstable. We discuss how variations in these parameters affect the evolution. We find that systems with equal-mass 
giant planets undergo the most violent instabilities, and that these destroy both terrestrial planets and the outer planetesimal disks that 
produce debris disks. In contrast, systems with low-mass giant planets efficiently produce both terrestrial planets and debris disks. 
A large fraction of systems with low-mass (M < 30 M e ) outermost giant planets have final planetary separations that, scaled to the 
planets' masses, are as large or larger than Uranus and Neptune in the Solar System. We find that the gaps between these planets are 
not only dynamically stable to test particles, but are frequently populated by planetesimals. The possibility of planetesimal belts be- 
tween outer giant planets should be taken into account when interpreting debris disk SEDs. In addition, the presence of ~ Earth-mass 
"seeds" in outer planetesimal disks causes the disks to radially spread to colder temperatures, and leads to a slow depletion of the 
outer planetesimal disk from the inside out. We argue that this may explain the very low frequency of > 1 Gyr-old solar-type stars with 
observed 24yum excesses. Our simulations do not sample the full range of plausible initial conditions for planetary systems. However, 
among the configurations explored, the best candidates for hosting terrestrial planets at ~ 1 AU are stars older than 0.1-1 Gyr with 
bright debris disks at 70/jm but with no currently-known giant planets. These systems combine evidence for the presence of ample 
rocky building blocks, with giant planet properties that are least likely to undergo destructive dynamical evolution. Thus, we predict 
two correlations that should be detected by upcoming surveys: an anti-correlation between debris disks and eccentric giant planets 
and a positive correlation between debris disks and terrestrial planets. 

Key words, planetary systems: formation — methods: n-body simulations — circumstellar matter — infrared stars — Kuiper belt 
— Solar System — astrobiology 



1. Introduction 

The Solar System's distinctive architecture, in which rocky ter- 
restrial planets lie interior to gas and ice giants, with the Kuiper 
Belt of smaller bodies beyond, is not unexpected. The mass 
in protoplanetary feeding zones increases with orbital distance, 
but the resulting tendency toward the formation of larger plan- 
ets further out is eventually frus trated both by the lengthen- 
ing ti me scale for accretion (e.g.. lLissaue n ll993HKokubo&Idal 
2002), and by the increased abili ty of planetary co res to scat- 
ter planetesimals inward dLevison & Steward |200 lb . The com- 
petition between these effects plausibly leads to the relatively 



slow (~100 Myr) assembly of a handful of terrestrial planets 
inside a few AU, the faster growth of several planetary cores 
with M > 5 M @ inside ~10 AU, and the persistence of a 
belt of unconsolidated debris further out. Simple arguments of 
this kind fail to establish how often planetary cores grow fast 
enough to admit the formation of fully fledged gas giants, but 
empirical estimates based on extrapolations of radial velocity 
and micr olensing surveys suggest that gas giant formation is 
common (ICumming et al . 2008; G ould et al.ll2010h . Debris disks 
dWvattll2008h are also observed around a significant fraction of 
young stars - despite the existence of both dynamical and col- 
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lisional processes that can destroy them on a time scale short 
compared to the main sequence lifetime of Solar-type stars - but 
the abundance of Earth-mass planets remains to be measured. 

For the Solar System, we have access to a unique array of 
observational constraints. Even with these advantages the ex- 
act nature of the interactions between the giant planets, the ter- 
restrial planets, and the Kuiper belt remain under debate. In 
the inner Solar System, at a minimum secular resonances with 
the g iant planets would have influenced terrestria l planet forma- 
tion ([Nagasaw a et al.1200 5: Raym ond et al.l2009d) . There could, 
however, have been stronger effects. Gas driven migration of the 
giant planets coul d have brought them (temporarily) closer to the 
Sun dWalsh et alJl201 ll) . directly reducing the supply of raw ma- 
terial i n the Mars re gion and preventing the growth of a larger 
planet dHansenll2009l) . In the outer Solar System, early studies 
focused on the dynamics of Neptune, which is of low enough 
mass that inward scattering of planetesimals can driv e substan- 
tial outward orbital migration dFernandez & Ipl l 1984). The mi- 
gration can deplete the mass in the Kuiper Belt and re sult in the 
reson ant capture of Pluto and other bodies by Neptune ( Ma lhotral 
1 1993b . Subsequent work introduced the idea of larger-scale 
dynamical instability a mong the Solar System's giant planets 
dThommes et al.l 1 1 999h . In the most-developed models, early 
outer Solar System evolution is characterized by a combination 
of planetesimal migration, close encounters be t ween planets, 
and r esonant interactions (Tsig anis et alj 120051 : ILevison et alj 
l20llb . 

Theoretically, attempts to construct equally detailed mod- 
els for extrasolar planetary system evolution are hampered 
by uncertainties in the distribution of the initial disk condi- 
tions, and by our po or knowledge of the evolution of gas 
disks ( Armitag g|201 lb and fo rmation mechanism for planetes- 
imals dChiang & Youdin 1201 Oh . Several observed properties of 
extrasolar planetary systems, however, including the existence 
of hot Jupiters (whose orbits are sometimes misaligned with re- 
spect to the stellar spin axis) and the prevalance of eccentric 
orbits, favor scenarios in which large-scale orbital evolution of 
giant planets is th e norm dWinn et alj|2010fc iTriaud et al.ll2010t 
Schlaufman 2010). It is therefore of interest to determine the di- 
versity of outcomes that could arise given initial conditions and 
dynamical processes similar to those of the Solar System, and 
to examine how those outcomes depend upon parameters such 
as the masses of the giant planets, and the properties of primor- 
dial planetesimal belts. Doing so is the goal of the present paper. 
We are particularly interested in studying how terrestrial plan- 
ets form, and debris disks evolve, in the presence of dynami- 
cally active giant pl anet systems. In an earlier paper ("Paper 1" 
iRavmond et al.ll20lTI) we showed that if giant planets form in or 
near dynamically unstable configurations, there are striking cor- 
relations between the nature of the terrestrial planets that form, 
and the properties of outer debris disks whose ongoing colli- 
sional evolution can be obs erved out to ages of several Gyr (e.g. 
Wvatt 20081 lKrivovll2010l) . The dynamically calm conditions 
that favor the formation of massive terrestrial planet systems 
also result in long-lived outer debris disks, that remain bright in 
cold dust emission (e.g. at A ~ lOfim) to late times. In systems 
that suffer more dramatic dynamical evolution, we identified a 
channel for the formation of unusual terrestrial planet systems 
in which a single planet exhibits large oscillations in eccentric- 
ity and inclination due to secular coupling to a scattered giant. 
Here, we consider a broader range of models within the same 
qualitative class, and study how robust our earlier conclusions 
are to changes in the poorly-constrained model parameters. 



Our paper is structured as follows. In section 2 an outline 
of our methods is presented (the reader is referred to Paper 1 
for more details and numerical tests). In subsequent sections we 
present results of different sets of simulations to test the effect 
of: the giant planet mass and mass distribution (section 3), the 
width and mass distribution of the outer planetesimal disk (sec- 
tion 4), and the presence of gas during giant planet instabilities 
(section 5). In section 6, we discuss the implications of our mod- 
els for debris disks and terrestrial planet systems. We conclude 
in section 7. 

2. Methods 

The initial conditions for our simulations assume that the loca- 
tion and mass of the rocky material in the terrestrial planet zone, 
and the mass of planetesimals in the outer disk, are fixed at val- 
ues similar to those employed in Solar System models. We as- 
sume that the masses of the giant planets, on the other hand, have 
a broad dispersion (extending up to masses above those realized 
in the Solar System) that is uncorrected with the mass in either 
the terrestrial planet region or in the outer planetesimal disk. All 
of our simulations contain three radially-segregated components 
orbiting a solar-mass star: 

1. The building blocks of terrestrial planets: 9 M ffi in 50 plane- 
tary embryos and 500 planetesimals from 0.5 to 4 AU, with 
equal mass in each component and a radial surface density 
profile S ~ r _1 . The initial eccentricities were chosen at ran- 
dom from 0-0.02 and the initial inclinations from - 0.5°. 

2. Three giant planets at Jupiter-Saturn distances: the inner- 
most planet is placed at 5.2 AU and the two others are 
spaced outward by 4-5 mutual Hill radii. We adopt three- 
planet initial conditions because this is the simplest plausi- 
ble configuration that evolves dynamically to match the mea- 
sur ed eccentricity distrib ution of massive extrasolar plan- 
ets dChatteriee et alj|2008l) . The planets were placed on ini- 
tially circular orbits with randomly-chosen inclinations of 
0-1°. 

3. An outer disk of planetesimals thought to be analogous to 
the primitive Kuiper belt. This belt consists of 1000 plan- 
etesimals with a total mass of 50 M ffi . The belt starts 4 Hill 
radii beyond the outermost giant planet and extends radially 
for 10 AU, and also follows an r _1 radial surface density pro- 
file, f The initial eccentricities were chosen at random from 
0-0.01 and the initial inclinations from - 0.5°. 

Adopting these initial conditions amounts to making im- 
plicit assumptions about the typical outcome of planet forma- 
tion. First, our terrestrial, giant planet, and outer disk zones 
are located such that they are in immediate dynamical contact 
with each other. This is reasonable only if planetesimal for- 
mation results in a smooth, gap-less distribution of bodies in 
0.5 AU < a < 20 AU, and if giant planet migration is lim- 
ited. Substanti al giant planet migra t ion, o f the kind envisaged 
in models by iMasset & Snellgrovel (1200 ll) . IWalsh et all (1201 ll) 
and lPierens & Raymond! d20 11 ), could create dynamical separa- 
tion between the giant and terrestrial planets prior to the gas-less 
phase of evolution that we simulate. Second, we assume non- 
resonant initial conditions for the giant planets. Mean-motion 
resonances can be established readily if there is significant mi- 
gration, due to either gas disk torques or planetesimal scatter- 
ing, and a plausible alternate class of models could be con- 
structed in which fully resonant initial conditions were the norm 
dMorbidelli et al.ll2.007b . We do not consider this possibility fur- 
ther here. 
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Embryo and giant planet particles feel the gravitational at- 
traction of all other bodies in the simulation. Planetesimal parti- 
cles, both in the inner and outer disk, feel the gravity of embryos 
and giant planets but do not self-gravitate. This commonly-used 
approximation allows for an adequate treatment of collective 
particle effects and dramatically reduces the required computa- 
tion time. Our methods are outlined in detail in Paper 1 . Here we 
summarize the key points. Each simulation was integrated for 
100-2 00 million years using the hybrid version of the Mercury 
code (Cham berslTT999h with a 6 day timestep. Collisions be- 
tween particles were treated as inelastic mergers. Particles were 
removed from a simulation when they either came within 0.2 AU 
of the central (solar-type) star at which point they are assumed 
to have collided with the star, or ventured farther than 100 AU 
from the central star (1000 AU for the widediskruns discussed 
in Section 4.2), at which point the particle is assumed to have 
been ejected from the system. 

In Paper 1 we presented the results from our fiducial mixed 
set of simulations. In these simulations the giant planet masses 
are d rawn randomly from the observed exop lanet mass distribu- 
tion dButler et alJl2006HUdrv & S antos 20071) : 



dN 
dM 



oc M 



(1) 



where masses were chosen between one Saturn mass and 
3 Jupiter masses. In these runs, the masses of individual plan- 
ets were chosen independently. 

We also computed a number of alternate models in which 
either the range of the mass function, the assumption of inde- 
pendent masses, or the properties of the outer disk were altered 
(see Table Q}: 

- The lowmass simulations represent systems with low-mass 
giant planets. For these cases the giant planet masses also 
follow the observed exoplanet distribution, but with masses 
between 10 M ffi and 1 Mj llp . The initial conditions for the gi- 
ant planets in these lowma ss simulations a r e the same as th e 
"mixe d2" simulations in iRavmond etall (l2008aL l2009blal 
2010). Abnormally low giant planet masses, relative to 
the mass in terrestrial planet-forming material, might occur 
physically in disks around stars where stronger than average 
photoevaporation limits the disk lifetime. 

- The equal simulations comprise four sets of simulations, 
each containing three giant planets with fixed masses of 
30M ffl , IMgat, lMj, or 3Mj. For the equal simulations the 
planets were placed in a slightly more compact configura- 
tion (separated by 3.5-4 mutual Hill radii rather than 4-5) 
to ensure that they would become unstable. This variation 
mimics the reality that the conditions that favor the growth 
of one gas giant to high masses - for example early core 
formation or a long disk lifetime - probably apply also to 
other planets in the same system. Past work suggests that 
thes e simulations should produce the most violent instabili- 
ties dFord et al.| [2003: Ra ymond et alj|2010l) . 

- The widedisk simulations test the effect of planetesimal 
disks that are 20 AU wide rather than 10, and twice as mas- 
sive (so that there is the same mass in the first 10 AU of the 
annulus as the fiducial case). In these simulations, the radius 
beyond which an object is considered ejected was 1000 AU 
rather than 100 AU. The giant planets' initial orbits are iden- 
tical to the mixed set. 

- The seeds simulations test the effect of the mass distribution 
within the planetesimal disk by including five or ten equally- 
spaced equal-mass fully self-gravitating seeds of either 2 M e 



or 0.5 M ffi respectively. The total mass and width of the plan- 
etesimal disk was held fixed at 10 AU. 
- The gas simulations test the effect of the presence of a 
gas disk during and after the giant planet instabilities. (The 
methodology is discussed in §5.) 

Each simulation was post-processed to calculate the spectral 
en ergy distribu t ion of dust in the system following the method 
of iBooth et ail d2009t) with a few small changes (see Section 
2.3 in Paper 1). To do this each planetesimal particle was as- 
sumed to represent a population of objects with sizes between 
2.2yum and 2000 km. This population was assumed to be in col- 
lisional equilibrium such that the differential size distribution 
can be written as n(D) oc D 3 5 dDohnanvil 1 1 9691) . The radial 
distribution of dust was calculated by a simple combination of 
the planetesimal orbital distribution (and by sampling eccen- 
tric orbits at multiple intervals along their orbit equally spaced 
in mean anomaly). The spectral energy distribution was calcu- 
lated by assuming that the dust grains in each radial bin emit as 
blackbodies based on their effective temperature. At each sim- 
ulation timestep the collisional timescale t c was calculated for 
the largest objects (D = 2000 km) for the population of both 
asteroidal and cometary planetesimals. In practice, t c represents 
the mean timescale between collisions energetic enough to dis- 
rupt an object of a given size at a given orbital radius. t c is a 
function of the mass, width, and orbital distribution of the plan- 
etesimal belt as well as the physical properties of planetesimals 
themselves, in particular their Q* D , the impact energy needed 
to catastrophically disrupt a planetesimal of size D (for details, 
seelWyatt et al.lll999ll20 07b; Boo th et alj|2009t IRavmond et al.1 
1201 lL iKains et aTl 1201 ll) . Once t c was calculated for a given 
timestep, the effective dust mass of each population was de- 
creased by a factor of [1 + t/t c (D c )] . This decrease in the dust 
mass is not self-consistent because the planetesimal mass in the 
simulations is constant. This effect can be important for the aster- 
oidal planetesimals because their collisional timescale, t c ~ 10 4 
years, is short compared to the interesting timescales for dy- 
namical evolution. The opposite ordering typically applies for 
the outer, cometary planetesimals, whose collisional timescale 
is t c > 10 8 years. The dust fluxes used in the analysis later in the 
paper are dominated by the cometary component so this incon- 
sistency has little to no effect on our results. In addition, for the 
case of low-mass giant planets that migrate due to planetesimal 
scattering (i.e., the lowmass simulations), the timescale for dy- 
namical mass loss from the outer planetesimal disk is roughly an 
order of magnitude shorter than the timescale for the calculated 
collisional mass loss in that same region. Thus, our assumption 
that the planet-planetesimal disk dynamics is not affected by the 
collisional cascade appears reasonable for the outer disk. 

This simple model is based on previous studies that fit the 
statistics of debris disks using models for the collisional evo- 
lution of planetesimals ([Dominik & Decinl 120031: iKrivoy et all 
| 2005l l2006t IWvatt et alJl2007bi IWvattll2008t lLohne et alJ l2008 ; 
IKains et all 1201 ll) . Our model agrees to within a factor of 2- 
3 at 24/im and 70^m with more detailed calculation of dust 
production d uring the collisional evolution (iKenvon & Bromlev 
2008, 2010) a nd also with dust fluxes observed around solar- 
type stars ( Habing et al. 2001; Beichm an et al.l20 06: Moor et al 
2006; [Trilling et alJl2008t iHillenbrand et alJl2008t iGaspar et al.. 
2009; Carpent er et al.l 20091) . However, due to our incomplete 



knowledge of the physical properties of planetesimals, there re- 
mains uncertainty in the dust fluxes of up to a n order of magni- 
tude for a given system (see Boo th et al.ll2.009b . Our model does 
not include outgassing from comets (i.e., outer disk planetesi- 
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Table 1. Summary of the simulations 



Set 


N(sims)* 


Giant planet mass distribution 


Giant planet spacing 


Gas? 


mixed 


156 


dN/dM ~ AT 11 from lM Sat to 3Mj up 


4 - 5Rh* 


no 


equal 


68 


M = 30M ffl , lM Sat , lM Jup or3M Jup 


2.5 - 2>R H , m 


no 


lowmass 


86 


dN/dM ~ M~ lA from 10M e to lM Jup 


Same as mixed 


no 


widedisk 


91 


Same as mixed 


Same as mixed 


no 


smallseed 


44 


Same as mixed 


Same as mixed 


no 


bigseed 


47 


Same as mixed 


Same as mixed 


no 


gas 


45 


Same as mixed 


crossing orbits 


yes 



*Note that N(sims) represents the number of simulations that met our criteria for run time of > 100 Myr and energy conservation of dE/E < 10~ 2 
(see Appendix A of Paper 1 for numerical tests). Our initial batches of simulations were somewhat larger: 200 simulations in mixed, 20 for each 
equal set, 100 for lowmass and widedisk, and 50 for gas, smallseed and bigseed. 



mals) when they enter the inner Solar System. Thus, the dust 
fluxes that we calculate during planetesimal bombardments are 
significantly underestimated. Indeed, the 24 ^im dust flux during 
the late heavy bombardment calculated by iBooth etail (2009) 
with our method reaches a peak tha t is roughly an order of mag- 
nitude lower than that calculated bv lNesvornv et alj (12010), who 
accounted for cometary dust production. Cometary outgassing 
is of importance at mid-infrared wavelengths (10 < A < 50/im) 
during and shortly after bombardments. As our results focus on 
the steady-state production of cold dust in outer planetesimal 
disks rather than on bombardments, we are not strongly affected 
by this effect. 

An important point is the fact that our sets of simulations sys- 
tematically over-predict the frequency of debris disks by a factor 
of roughly 2-4. Given Spitzer's detection limits, the observed 
frequency of debris disks at 70/zm around Solar-type stars olde r 
than 1 Gyris 16.4% dTrilling et alj|2008t ICarpenter et al.ll2009l) . 
At 24yum, the observed freq uency is much lower, only 2-3% for 
stars older than 300 Myr (ICarpenter et al.1 l2009t iGaspar et all 
2009). Various aspects of and potential solutions to this issue 
will be discussed in more detail throughout the paper. 

2.1. An example mixed simulation 

Here we briefly present a simulation from the mixed set to fa- 
cilitate comparison against example simulations from other sets 
that are presented in later sections. The chosen system started 
with giant planets of 0.96, 0.46 and 0.64 Mj in order of increas- 
ing orbital distance, and an outer planetesimal disk that extended 
out to 22.8 AU. 

Figure [TJ shows the evolution of the simulation's dynamics 
and calculated dust flux. The system remained stable for 42.8 
Myr at which time it underwent a strong dynamical instability 
that started with a close encounter between the middle and outer 
giant planets that triggered a series of planet-planet scattering 
events over the next 400,000 years. The instability culminated 
in the ejection of the middle giant planet, and the surviving two 
giant planets swapped orbits (i.e., the innermost planet became 
the outermost and vice versa). At the end of the simulation both 
planets' eccentricities are large: e,, me . r oscillates between 0.61 
and 0.83 and e outer between 0.15 and 0.37. Despite their large ec- 
centricities, both of the planets' inclinations with respect to the 
initial orbital plane remain modest (at least in comparison with 
the planets' eccentricities): im mr oscillates between 3.5° and 14° 
and iouter between zero and 3.1°. With a semimajor axis of 2.55 
AU, the orbit of the inner giant planet is well-represented by the 
more eccentric of the known exoplanets, while the outer planet 
would probably not be currently detectable. 

Before the instability the inner planetary system was under- 
going standard terrestrial accretion. Embryos' eccentricities re- 



mained damped by dynamical friction from planetesimals and 
embryos grew by frequent planetesimal impacts and occasional 
giant embryo impacts. At the time of the instability the accretion 
process was relatively mature, as only eight embryos remained 
inside 2 AU with masses between 0.08 and 0.91 M ffi , as well as 
three smaller (0.07 - 0.2 M ffi ) embryos in the asteroid belt. In the 
immediate aftermath of the instability, ten of the eleven embryos 
collided with the central star. The mechanism that drove the em- 
bryos into the star was strong eccentricity pumping by a combi- 
nation of close encounters and secular pumping by the (initially 
outermost) inward-scattered giant planet. The one embryo that 
was not driven into the star was the outermost one, which had a 
pre-instability semimajor axis of 2.6 AU and was ejected after a 
series of close encounters with the scattered giant planets. 

The pre-instability outer planetesimal disk was for the most 
part dynamically calm. The inner edge of the disk was slowly 
eroded during this period as the giant planets cleared out plan- 
etesimals that were unstable on long timescales. A few other 
resonances (e.g., the 3:1 resonance with the outer giant planet 
at 21.4 AU) also acted to increase the eccentricities of long-term 
stable planetesimals. A slow trickle of planetesimals was also 
destabilized by certain strong mean motion resonances, notably 
the 2:1 resonance. When the instability started, the (previously 
middle) giant planet was scattered out into the planetesimal disk 
on an eccentric orbit. Its eccentricity was further pumped by a se- 
ries of close encounters with the (initially innermost) planet such 
that for a period of several hundred thousand years (several thou- 
sand cometary orbits) the giant planet's orbit completely crossed 
the initial planetesimal disk. The outer planetesimal disk was 
entirely destabilized by secular interactions and close encoun- 
ters. These planetesimals' eccentricities increased drastically un- 
til they either hit the sun (this occurred about 25% of the time), 
were scattered out beyond 100 AU and removed from the sim- 
ulation by presumed hyperbolic ejection (~ 75%), or collided 
with a giant planet (~ 0.5%). More than 80% of the outer disk 
planetesimals were destroyed within 500,000 years and 97% 
within 5 Myr. Only 23 planetesimal particles survived more than 
5 Myr after the instability on orbits that were unstable on 10 
Myr timescales, typically with high inclinations and eccentrici- 
ties. At the end of the simulation a single planetesimal survived, 
although it is almost certainly unstable on longer timescales as 
its orbit crosses the outer giant planet's. 

The system's spectral energy distribution (SED) - shown at 
right in Figure [TJ - reflects its dynamical evolution. The aster- 
oidal planetesimals are quickly ground to dust, as their colli- 
sional timescales are only 10 4 - 10 5 years. This causes a rapid 
decline in flux at short wavelengths (A < 20yum). The erosion of 
the inner edge of the outer planetesimal disk causes a continued 
decrease in flux at shorter wavelengths. The flux at A > 50/im 
is dominated by the mass in outer disk planetesimals and is only 
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Fig. 1. Evolution of a reference mixed simulation in which the giant planets underwent a violent instability after 42.8 Myr of evolu- 
tion. The initial giant planet masses were, in order of increasing orbital distance, 0.96, 0.46 and 0.64My. Left: Orbital eccentricity 
vs. semimajor axis of each body in the simulatio n. The size scales w ith the mass 1/3 and the color corresponds to the water content, 
using initial values taken to Solar System data dRavmond et al.ll2004l) and re-calculated during impacts by mass balance. The giant 
planets are shown as the large black bodies and are not on the same size scale. Top right: The spectral energy distribution of the dust 
during five simulation snapshots. The dashed line represents the stellar photosphere. Bottom right: The ratio of the dust-to-stella r 
flux F jF star at 25yum as a function of time. The rough Spitzer observational limit is shown with the dashed line (Trilli ng" et al.l l2008). 
A movie of this simulation is available at http://www.obs.u-bordeauxl.fr/e3arths/raymond/scatterSED_199.mpg 



very weakly affected by the grinding of asteroids or the slow 
erosion of the inner edge of the planetesimal disk. When the in- 
stability occurs and the outer planetesimal disk is destabilized, 
a large number of planetesimals are temporarily placed on high- 
eccentricity (and therefore small-periastron) orbits which pro- 
duce hot dust. This burst in hot dust changes the shape of the 
SED by increasing the flux at short wavelengths. This is the 
cause of the spike in flux seen at 25jum in the lower right panel of 
Fig-CD As noted already, the magnitude of the spike is underesti- 
mated because we do not acco unt for the outgassing that occurs 
when an icy body is heated (see lNesvornv et 311201 0*). However, 
as the outer planetesimal disk is removed the flux drops dramat- 
ically. The dust flux in the 30 Myr after the instability is main- 
tained at a relatively high level by a single planetesimal that sur- 
vived from 43 to 66.5 Myr between the giant planets with a per- 
ihelion that dropped periodically below 3 AU (on a retrograde 
orbit). This single close-in planetesimal produced enough dust 
to keep the system above the 25/im detection threshold during 
this period. We note that our dust production scheme does not 



allow for the collisional evolution of a single planetesimal so its 
dust flux is certainly overestimated. 

This example is relatively extreme in terms of the plan- 
ets' final orbital eccentricities and in that all the terrestrial and 
cometary particles were destroyed. However, as we will see be- 
low, this simulation allows for a convenient comparison with up- 
coming examples because the instability is delayed and so the 
evolution of the dust flux from the quiescent outer disk is unper- 
turbed at early times. 

3. Effect of the giant planet masses and mass 
distribution 

We now analyse two sets of simulations that explore alternate 
giant planet mass distributions than the mixed set analysed in 
Paper 1. The mixed set included three planet systems with the 
masses of the planets being chosen randomly and independently 
in the range between a Saturn mass and 3 Jupiter masses. In the 
equal simulations, planet masses within a given system are the 



5 



Raymond et al.: Debris disks as signposts of terrestrial planet formation II 



mixed lowmass 




'. i i i i i i i i i i i : c i i i i i i i i i i i ; 

2 3 5 10 20 30 50 2 3 5 10 20 30 50 

Semimajor Axis (AU) Semimajor Axis (AU) 

Fig. 2. Eccentricity vs. semimajor axis for the surviving giant planets in the mixed (left panel) and lowmass (right) simulations. 
Black circles represent unstable simulations - defined as simulations in which at least one giant planet-planet scattering event 
occurred - underwent and grey dots are stable simulations. The size of each circle is proportional to the logarithm of the planet 
mass. The dashed vertical lines represent the median outer edge of the planetesimal disk for each set of simulations. Note that the 
outer edge varied from simulation to simulation depending on the giant planet masses, from 20.6 to 27 AU (with a median of 23.7 
AU) for the mixed simulations and 17.3 to 23.3 AU for the lowmass simulations (median of 19.7 AU). 



same. We consider masses between 30 M© and 3Mj for different 
systems. Assuming that the masses within individual systems 
are perfectly correlated has the effect of maximizing the strength 
of dynamical instabilities, compared to systems where there is a 
range in masses. In the lowmass simulations, planet masses are 
drawn randomly from the observed distribution but only in the 
range of 10M© to IMj (in contrast to the range of Ms at to 3Mj 
for the mixed set). The inclusion of the lower mass planets in- 
creases the fraction of systems for which dynamical interactions 
between planets and the planetesimal disk are important. These 
interactions can take two forms. First, "planetesimal-driven mi- 
gration" changes the orbital radius of a planet due to the back- 
reaction of planetesimals that are gravitationally scattered by the 
planet, and thus changes a planet's orbital radius while maintain 
ing a small eccentricity 
Hahn & Malhotra 199S 



((Fernandez & Idll984t iMalhotr 



Gom es et alj l2004t iKirsh et al 



1993 



2009; 



Levison et alJ l2010). Second, the orbit of an eccentric planet can 
be re-circularized by "secular friction", a process by which an 
eccentric planet excites the eccentricities of the outer disk plan- 
etesimals an d causes a corresponding decrease in the planet's 
eccentricity (Thommes et al. 19991 iLevison et al 200$. A low- 
mass planet can therefore be gravitationally scattered by an- 
other planet in the inner part of a planetary system and have 
its eccentricity decreased on a much wide r orbit by secular fric- 
tion with the outer pla netesimal disk (e.g., Thommes et al. 1999; 
iRavmond et ai1l2010l) . 



The low mass and equal sets of simulations are of particu- 
lar interest because a combination of the two can produce an al- 
ternate sample that matches the observed exoplanet distribution. 
We explore the two sets of simulations independently (sections 
3.1 and 3.2) and later combine them into a sample to compare 
with exoplanet statistics (called case B in section 6). 



3.1. Systems with low-mass giant planets (the lowmass 
simulations) 

Figure [2] shows that, because of planetesimal-driven migration 
and secular friction, surviving low-mass giant planets populate 
different regions of parameter space than more massive giant 
planets. Given that all of our simulated giant planets started in 
the 5-15 AU region, massive giant planets are only able to al- 
ter their semimajor axes by energy exchange during close en- 
counters to essentially follow curves with perihelia or aphelia 
at the encounter distance. In other words, high-mass planets at 
large a necessarily have large e. However, low-mass planets can 
have large a and small e by either 1) planetesimal-driven migra- 
tion, which maintains planets' low e out to large a (to the outer 
edge of the planetesimal disk), or 2) being scattered outward in 
a dynamical instability but having their eccentricities damped 
by secular friction with the outer planetesimal disk. Similarly, 
massive planets that are scattered interior to 5.2 AU necessar- 
ily have large e but low-mass inner giant planets can undergo 
inward changes in a and survive on low-e orbits. Indeed many 
low-mass planets do just that, ending up at a = 2.5 - 4 AU. 
Note that secular friction is only relevant in unstable systems in 
which a planet acquires a large orbital eccentricity. On the other 
hand, planetesimal-driven migration is mainly relevant for sta- 
ble systems although periods of migration may occur in some 
cases after secular friction has already re-circularized the orbit 
of a scattered low-mass giant planet. 

Thus, the surviving high-mass giant planets retain a memory 
of their initial conditions: the stable planets and many unstable 
planets are clustered at their original locations. However, given 
the ease and inevitability of planetesimal-driven migration and 
secular friction for low-mass giant planets, the initial conditions 
are erased. 

Figure [3] shows the evolution of a lowmass simulation with 
initial giant planet masses of 12.4 M© (inner), 18.6 M© (mid- 
dle), and 35.9 M© (outer). In this simulation the giant planets 
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Fig. 3. Evolution of a simulation with low-mass giant planets. Left: Orbital eccentricity vs. semimajor axis of each body 
in the simulation. The size scales with the mas s 1 / 3 and the color corresponds to the water content, using initial values 
taken to Solar System data dRavmond et al.l 120041) and re-calculated during impacts by mass balance. The giant planets are 
shown as the large black bodies and are not on the same size scale. Top right: The spectral energy distribution of the 
dust during five simulation snapshots. The dashed line represents the stellar photosphere. Bottom right: The ratio of the 
dust-to-stellar flux at 25 microns as a function of time. T he rough observatio nal limit of the MIPS instrument on NASA's 
Spitzer Space Telescope is shown with the dashed line (Trilli ng st alJ [2008). A movie of this simulation is available at 
http://www.obs.u-bordeauxl.fr/e3arths/raymond/scatterJowmass_54.mpg 



underwent an instability after 33,000 years, which threw the 
inner planet into the outer planetesimal disk. Once interacting 
with the planetesimal disk, the planet scattered planetesimals in- 
ward and migrated outward for roughly 20 Myr, then slowed 
when it came within ~3.5 Hill radii of the outer edge of the 
planetesimal disk. As th is represents the approximate boundary 
for dynamical stability dMarchal & Bozislll982_[Gladmanll 19931: 
Chambers et al.|[l996l) . the number of planetesimals available to 
be scattered by the planet decreased and the planet's migration 
slowed drastically. The planetesimals that were scattered inward 
from the outer disk were for the most part subsequently scattered 
outward by the inner giant planets, causing the two inner plan- 
ets to migrate inward. However, some of the inward-scattered 
planetesimals were trapped on low-eccentricity orbits between 
the two inner giant planets and the outer giant as the outer planet 
migrated outward by continued planetesimal scattering. This is 
similar to the mechanism that may have been responsible for 
populating the Solar System's ast eroid belt during the outward 
migration of Jupiter and Saturn (Wals h et al.ll201 lh . The inner 



giant planets' inward migration was fueled by the planetesimals 
that ended up with high-eccentricity, low-perihelion orbits after 
being scattered inward, although the inner giants also scattered 
some embryos and planetesimals from the inner disk. The inner 
giant planet migrated in to 2.69 AU but maintained an eccentric- 
ity lower than 0.1 throughout and less than 0.05 during the last 
phases. 

Two terrestrial planets formed in the simulation shown in 
Figure [3] a 1.25 M ffi planet at 0.61 AU and a 0.69 M e planet 
at 1.11 AU. The eccentricities of the inner and outer planet are 
0.06 and 0.1 1, respectively, with peak to peak oscillation ampli- 
tudes of 0. 1 1 and 0.20. The inner planet underwent its last giant 
(embryo) impact after 40.6 Myr but the outer planet did not un- 
dergo any giant impacts after 3.8 Myr. The inner planet is wet: 
it accreted a small amount of water from material that originated 
in the inner asteroid belt (an embryo and two planetesimals from 
~ 2 AU) but the bulk of its water came from a single cometary 
impact. The outer planet did not accrete any m aterial from be- 
yond 2 AU and so is considered to be dry (see Raymond et al.l 
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2004, 120071) . The closest giant planet's semimajor axis is only 
1.6 AU larger than the outer terrestrial planet's but given the gi- 
ant planet's small mass (19.2 M ffi ) and modest eccentricity (0.035 
with oscillations of 0.02 in full amplitude) the system is stable. 

In the simulation from Fig. [3] only a relatively small frac- 
tion of the initial terrestrial mass was incorporated into the two 
surviving terrestrial planets. The majority of the initial terrestrial 
mass (57%) was ejected from the system and an additional 10% 
collided with the central star. This contrasts with the unstable 
systems with higher-mass giant planets, in which terrestrial ma- 
terial preferentially collides with the star (as in the simulation 
from Fig. [TJ. This difference is due to the fact that more mas- 
sive giant planets pump the eccentricities of terrestrial bodies 
efficiently and can thus drive down their perihelion distances on 
short timescales. Low-mass giant planets require longer to excite 
high eccentricities and also migrate in reaction to scattering such 
bodies, inward in this scenario. Thus, in systems with unstable 
low-mass giant planets terrestrial material is more easily trans- 
ported outward than inward, to be ejected after many encounters 
with the outer giant planets. 

In addition, in the simulation from Fig.[3]two Mars-sized em- 
bryos were scattered out and survived on distant orbits, at 29.2 
and 34.8 AU, and one embryo collided with the middle giant 
planet. 

At the end of the simulation there are two surviving plan- 
etesimal belts: a low-eccentricity belt between the two inner and 
the outer giant planets and an outer disk of higher-eccentricity 
objects exterior to the outer giant planet. The outer belt is anal- 
ogous to the Solar Sys tem's scattered disk dLuu et alJ Il997t 
Dunc an & Levisonl 19971) . having been scattered by the outward- 
migrating giant planet. The scattered belt contains 3.7 M© in 72 
particles with a median eccentricity of 0.27 and a median incli- 
nation of 11.5°. This scattered disk also contains two embryos 
from the inner disk, including one that originated inside 2 AU. 
The inner belt of planetesimals - located between roughly 8 and 
14 AU - contains 1.3 M© in 27 planetesimals with a median ec- 
centricity of 0.12 and a median inclination of 16.8°. The orbital 
distributions and surface densities of these two populations are 
quite different, and we suspect that a wide diversity of planetes- 
imal belt structures must exist around other stars. 

The evolution of the system's dust brightness is shown in 
Figure[3] The SED of the system decreases systematically as the 
system loses mass, but changes shape after roughly 80 Myr when 
four separate icy planetesimals entered the very inner planetary 
system and remained on orbits interior to the innermost giant 
planet (with perihelion distances as small as 0.3 AU) for several 
tens of Myr before being ejectedQ At wavelengths longer than 
~ 50^im, the dust brightness decreased monotonically in time. 
However, shorter wavelengths (such as 25/itn; Fig. [3]) show the 
additional structure caused by the icy planetesimals entering the 
inner planetary system because they are sensitive to hot dust. 

As a whole, the lowmass simulations were extremely ef- 
ficient at forming terrestrial planets and also at creating long- 
lasting debris disks. Out of the 86 total simulations, 82 (95.3%) 
formed terrestrial planet systems containing a total of at least 
0.5 M e . Of the four remaining systems, three destroyed their 
terrestrial planets entirely and the fourth formed a single planet 



lowmass 




1 Note that the plateau in brightness seen in the 25yum plot of Fig. [3] 
may be slightly overestimated because our method does not adequately 
account for collisional grinding of bodies on extremely close-in orbits 
when they are isolated particles. In this case, these four isolated parti- 
cles were at close enough distances to dominate the flux at wavelengths 
shorter than ~ 50/um from about 80-120 Myr. 



1 2 3 4+ 
Number of terrestrial planets (>0.05 M Earth ) 

Fig. 4. Distribution of the number of surviving terrestrial planets 
for the unstable (grey) and stable (dashed) simulations in the 
lowmass set of simulations. 



of 0.25 M e . Of the 86 lowmass simulations, 76 (88.4%) were 
above the Spitzer detection threshold at 70/mi after 1 Gyr (73 
remained bright after 3 Gyr; recall that this is far higher than 
the observed frequency of 16.4%). Of the 54 (62.7%) unstable 
simulations, only 4 did not yield at least 0.5 M© in terrestrial 
planets and only 10 were not detectable at 70/im after 1 Gyr. 
The few systems with destructive instabilities were those that 
by chance contained several massive planets, and so essentially 
overlapped with the mixed distribution. Figure [4] shows that all 
32 stable simulations finished with two or more terrestrial plan- 
ets and also with bright debris disks. Note that in this figure we 
use a low mass cutoff of just 0.05 M© in our definition of a ter- 
restrial planet. Any surviving planetary embryo is therefore con- 
sidered a planet. This allows for a consistent comparison with 
other sets of simulations including more violent instabilities (like 
the equal simulations) in which surviving embryos are com- 
mon. In the Solar System, Mars is thought to be a surviving em- 
bryo dDauphas & Pourma ndl201 lb . 

Among systems with planets less massive than 50-100 M© 
(~ 0.5 - I Ms at), there was little difference in the final out- 
come between systems that underwent planet-planet scattering 
and those that did not. This is because the planetesimal disk pro- 
vides strong enough damping to quickly decrease the planets' 
eccentricities back to near zero. Systems containing a single rel- 
atively massive giant planet (M > Ms at ) also ended in a dynami- 
cally calm state because instabilities caused the lower-mass giant 
planets to be scattered and, again, their eccentricities and incli- 
nations are quickly damped. The only situation that preserved 
large eccentricities was the relatively infrequent combination of 
multiple massive giant planets in the same system. In those sys- 
tems the large eccentricity caused by strong scattering between 
giant planets could not be damped (the low-mass giant planets 
in such systems are usually scattered, sometimes to be ejected or 
sometimes re-circularized in the outer planetesimal disk). Thus, 
unlike massive giant planets, the eccentricities of low-mass plan- 
ets do not retain a memory of the system's dynamical history. 
In addition, multiple giant planets must exist in the same sys- 
tem to yield eccentric g iant planets. The abundance of observed 
eccentric planets (e.g. lWright et al.l[2009t) thus points to the fre- 



8 



Raymond et al.: Debris disks as signposts of terrestrial planet formation II 




Mass of outermost giant planet (Earths) Width of stable zone A (mutual Hill radii) 

Fig. 5. In the lowmass planetary systems, the fraction of stable zones between pairs of outer giant planets that contain at least one 
planetesimal on a stable orbit as a function of the mass of the outermost giant planet (left panel) and the width of the stable zone A in 
units of mutual Hill radii (right panel). For the left panel, bins were evenly logarithmically spaced from 10 M ffl to IM j. For the right 
panel, bins were evenly spaced from A = 3.5 to 30. The error bars were calculated using binomial statistics followingFBurgasse r et al] 
(2003). The dashed line marks the two planet stability boundary (iMarchal & Bozisll 19821 iGladmanll 19931) . 



qu ency of strong instabilities (although alternate models exist; 
see Ford & Rasio 2008, for a thorough review). 

Most of the lowmass systems have large gaps between gi- 
ant planets, and many of these gaps contain planetesimals on 
stable orbits, such as the belt between 8-14 AU in the system 
from Fig. [3] The existence of isolated belts alters the radial dis- 
tribution of dust and can be inferred from the spectral energy 
distribution. There are several systems known to contain belts of 
dust that appear to be radially confined, presumably by known or 
as-yet undetected pla nets (e.g., Beich man et al.l20 05: Lisse et all 
l200aiSuetall2009T) . 

The majority of the lowmass systems (60/85 = 71%) have 
gaps between the outer two giant planets with separations A of 
10 or more mutual Hill radii (in one system just a single gi- 
ant planet survived and so is not counted). Given that two plan- 
ets must be separated by at least A > 2"y3 » 3.5 mutual Hill 
radii for long-term dynamical stability (Marchal & Bozis 1982; 
lGladman|[l99l . the existence of a zone between two planets' 
orbits that is stable for planetesimals requires at a minimum an 
interplanetary separation A > 10. In practice, somewhat wider 
gaps are needed in realistic systems. The separation between 
Saturn and Uranus, and between Uranus and Neptune, amounts 
to 14 mutual Hill radii, and there is only a small region that 
is stable over long timesc ales between the latter two planets 
dHolman & Wi sdom 1993). In our runs we frequently find gaps 
that are not only wider than their Solar System counterparts, but 
which are also populated with primordial material despite the 
relatively coarse sampling of the outer planetesimal disk popula- 
tion. The lowmass simulations yield a range of separations from 
A = 3.9 to 30 and outer giant planet masses of 1 1 M e to 0.97 Mj. 
From the entire sample, the stable zone between the outer two gi- 
ant planets contained at least one particle in half of all lowmass 
simulations (38/85 = 45%, although 5 stable zones contained 
just a single planetesimal). The closest separation between two 
planets for which planetesimals existed on stable orbits between 
the two was A = 11.9, and the widest separation for which no 
planetesimals existed between two planets was A = 21.5. The 



total mass in these planetesimal belts ranged from the mass of 
one planetesimal particle (0.05 M ffl ) to 5.9M®, and in one case 
an embryo from the inner disk survived in such a belt. 

The two most important factors that determine whether a sta- 
ble zone between outer giant planets contains a planetesimal belt 
are first, the width of the stable zone and second, the mass of the 
outermost giant planet. Figure [5] shows that the probability that 
a stable zone contains at least one planetesimal particle on an 
orbit that is stable for long timescales ( i.e., whose orbit does not 
come within 4 Hill radii of any giant planet's orbit) as a func- 
tion of the mass of the outermost giant planet and the width 
of the stable zones (as quantified by A). Stable zones are pref- 
erentially filled for larger A values simply because there is a 
larger region of parameter space into which planetesimals can 
be scattered and survive, and the fraction of stable zones that is 
filled increases dramatically for A > 15. Stable zones are also 
preferentially filled in systems with outermost giant planets less 
massive than ~ 50 M e , and almost 100% of stable zones are 
filled when the outermost giant planet is less than about 20 M ffl 
(Flg-EJl. As they interact with the outer planetesimal disk, lower- 
mass planets scatter planetesimals onto lower eccentricities than 
do higher-mass planets, and these planetesimals are more likely 
to avoid encountering giant planets and to remain on stable orbits 
than if their obits are more eccentric. 

Do other system parameters influence the probability that a 
stable zone between giant planets will contain planetesimals? 
We tested the importance of several other parameters using a 
suite of Kolmogorov-Smirnov (K-S) tests that compared differ- 
ent characteristics of stable zones with and without planetesi- 
mals. Using a cutoff of p < 0.01 for a statistically significant dif- 
ference between the two populations, we found that five param- 
eters influence whether a stable zone will contain planetesimals. 
In order of most important (lowest p value) to least important 
(highest p value), these are 1) the width of the stable zone (A), 
2) the mass of the outermost giant planet, 3) the semimajor axis 
of the outermost giant planet, 4) the mass ratio of the outer two 
giant planets (systems with a large inner/outer mass ratio pref- 
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Fig. 6. Distribution of the number of surviving terrestrial planets 
for the equal set of simulations. Systems with giant planets of 
Msm or larger are shown in grey, and systems with giant planets 
of 30 M ffl are shown by the black dashed line. 



erentially contain planetesimal belts), and 5) the total mass in 
the surviving giant planets (the probability of stable zones being 
empty increases for higher total mass). The effect of these five 
parameters is statistically significant, although several are corre- 
lated (e.g., A correlates with the semimajor axis of the outermost 
giant planet, although A is about three orders of magnitude more 
important in determining whether a stable zone will be filled). 
We tested five additional parameters whose effect turned out to 
be unimportant (p > 0.1 in each case): the mass and semimajor 
axis of the giant planet marking the inner boundary of the sta- 
ble zone, the eccentricity of the outermost giant planet, and the 
mass-weighted eccentricities of all surviving giant planets. 

3.2. Systems with equal-mass giant planets (the equal 
simulations) 

While the lowmass simulations represent a calm environment 
conducive to the production of both terrestrial planets and 
bright debris disks, the equal simulations were destructive on 
both counts. This is simply because scattering among equal- 
mass giant planets is th e most violent planet-planet instabil- 
ity dRavmond et al.ll2~010l) . and strong instabilities destroy small 
bodies in both the inner and outer disks, typically by driving 
a large fraction of inner bodies into the central star and eject- 
ing the majority of outer bodies. Indeed, the eccentricity dis- 
tribution of the surviving equal giant planets with masses of 
Msot or larger was skewed toward higher values than the mixed 
sample, reaching values as high as 0.89 and with a median of 
0.35 (compared with 0.21 for the mixed giant planets). An ex- 
ception to this rule are systems with equal-mass giant planets 
that are themselves low-mass. In systems containing three gi- 
ant planets of 30 M®, the outcome was similar to the lowmass 
systems because planet-disk interactions trumped planet-planet 
scattering. Indeed, the low-mass equal systems were dominated 
by planetesimal-driven migration of the outer giant planet and 
behaved very similarly to the lowmass simulations. Given this 
strong dichotomy, we now consider just the high-mass equal 
simulations, as the lower-mass cases are more appropriately in- 
cluded with the lowmass set. 



Figure [6] shows that, of the 53 simulations with giant planet 
masses of Ms„t, Mj, or 3Mj, 36 (68%) destroyed all terres- 
trial material. Five simulations (9%) formed a single terres- 
trial planet, although these were all small, roughly Mars-mass 
~ 0.1 M e planets, and 4/5 of these were lone surviving em- 
bryos. The remaining 12 simulations with high-mass giant plan- 
ets (21%) formed two or more terrestrial planets. In contrast, not 
a single simulations with a 30 M e giant planet destroyed all of its 
terrestrial material. Rather, these systems usually formed several 
terrestrial planets - only 3/15 systems formed just one. 

Of the 53 equal systems, only 17 (32%) contained de- 
tectable amounts of cold dust at 70yum after 1 Gyr (only 14/53 
= 26% after 3 Gyr, although again note that this is higher than 
the observed frequency of 16.4%). Of the 17 systems with de- 
tectable dust at 1 Gyr, 11 (65%) formed terrestrial planets. Of 
the 36 systems with no detectable dust, only 6(17%) formed ter- 
restrial planets. Thus, there exists a natural connection between 
debris disks and terrestrial planets that spans the domain of gi- 
ant planet mass. This same debris disk-terrestrial planet correla- 
tion was found in Paper 1 for the mixed simulations. The details 
of this correlation depend on the system parameters and are of 
course confined to the context of our initial conditions, in partic- 
ular to systems with relatively massive outer planetesimal disks 
(see discussion in Section 6). 

The equal systems are violently unstable by construc- 
tion and so the outcomes tend to be extreme. This is 
precisely the type of behavior that is required by exo- 
planet observations, in particular the trend for more mas- 
sive planets to have more eccentric orbits than lower- 
mass planets d Jones et al . 2006; R ibas & M iralda-Escude 2007; 
Wrigh t et all 120091: iRavmond et alJ l2010h . If planet masses 
within individual planetary systems were random (as in the 
mixed simulations), then lower-mass giant plan ets should have 
highe r eccentricities than higher-mass giants ( Raymon d et al.1 
2010); this is the opposite of what is observed. Thus, the equal 
simulations provide a key ingredient in constructing a sample of 
simulations that matches the observed giant exoplanets, as dis- 
cussed in Section 6. 



4. Effect of the properties of outer planetesimal 
disks 

We now turn our attention to the effect of the properties of outer 
planetesimal disks. We first examine the seeds simulations - 
subdivided into the bigseed and smallseed sets - that con- 
tained a population of ~Earth-mass embryos in their outer plan- 
etesimal disks. We then test the effects of doubling the width 
(and total mass) of the planetesimal disk in the widedisk simu- 
lations. 



4. 1. The mass distribution of the planetesimal disk (the 
seeds simulations) 

In the seeds simulations a small number of fully-interacting 
massive bodies were included in the outer planetesimal disk 
(in contrast to planetesimal particles, which interact gravitation- 
ally with massive bodies but not with each other). The disk 
maintained the same total mass and numerical resolution (the 
masses of individual planetesimals were decreased to main- 
tain a constant total mass). In 50 simulations comprising the 
bigseed set, five icy embryos of 2M e each were included at 
11, 13, 15, 17 and 19 AU. In 50 additional smallseed simula- 
tions 10 seeds of 0.5 M e each were spaced with 1 AU of sep- 
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aration from 10.5 to 19.5 AU. These seed masses are broadly 
consis tent with calculations of accreti on in outer planetesimal 
disks dKenvon & Bromlevll2008ll2010l) . 

Figure|7]shows the evolution of one bigseed simulation that 
became unstable after a long delay and therefore allows for a 
comparison with the evolution of both stable (prior to the insta- 
bility) and unstable systems. The evolution is qualitatively sim- 
ilar to the mixed simulations in the inner disk such as the sys- 
tem from Fig.Q~| as embryos maintain smaller eccentricities than 
planetesimals via dynamical friction and accrete as their feed- 
ing zones begin to overlap. The instability started after 55.35 
Myr with a close encounter between the two outer giant plan- 
ets. After a series of close encounters lasting 40,000 years, the 
middle giant planet was ejected. At the end of the simulation, 
the system contains two giant planets each with e * 0.3 and 
two roughly Earth-sized terrestrial terrestrial planets (details in 
caption of Fig. |7]l. 

The evolution of the outer disk differs significantly from sim- 
ulations without seeds. The presence of the seeds introduces an 
effective viscosity into the outer planetesimal disk that is driven 
by scattering of planetesimals by close encounters with the icy 
embryos. This causes the disk to spread out radially. Within 1 
Myr the outer edge of the low-eccentricity portion of the disk 
has expanded from 22 to 26 AU, and beyond 30 AU in the next 
few Myr (Pig.[7J. This outward expansion is balanced by the loss 
of a large portion of the total disk mass that is scattered inward to 
encounter the giant planets and be ejected from the system. The 
planetesimal disk is further depleted on a 20-50 Myr timescale 
by instabilities in the system of icy embryos that increases their 
eccentricities - these instabilities are analogous to but far weaker 
than instabilities between the giant planets (and are weaker still 
in the smallseeds systems). When the giant planets go unsta- 
ble the bulk of the planetesimal disk is rapidly ejected and the 
last icy planetesimal is removed just after 100 Myr. 

The dust production in the seeds simulations is also dif- 
ferent than in simulations with calmer planetesimal disks. The 
planetesimal population closest to the giant planets - the inner 
several AU of the outer planetesimal disk - is depleted in a few 
Myr as the disk "viscously" spreads. Thus, warm dust has a far 
shorter lifetime than in simulations with no seeds. This is re- 
flected in the evolution of the 25pm flux in Fig. [7] that drops 
below the detection threshold after ~ 45 Myr while the system 
is still stable. In contrast, in the example mixed (Fig. [TJ and 
lowmass (Fig. [3} simulations, the 25/itn flux was more than an 
order of magnitude higher after 40-50 Myr of evolution. 

The shortened lifetime of warm dust is reflected in its chang- 
ing spectral energy distribution (Fig. |7). Compared with simu- 
lations without seeds, there is a much faster decrease in flux at 
A < 100/im at early times as the region just exterior to the outer- 
most giant planet is cleared much more efficiently and to a larger 
radial separation. At longer wavelengths the flux also decreases 
more rapidly due to depletion by scattering from icy embryos, 
although at long wavelengths this is counteracted in part by the 
expansion of the dust disk to larger radii and therefore increased 
surface area (though lower temperature). There are spikes in the 
flux (seen at 25/itn) when objects enter the inner planetary sys- 
tem to encounter and be ejected by the innermost giant planet. 
These spikes are due to the fact that each particle represents a 
distribution of smaller bodies and, with a higher numerical reso- 
lution, these spikes would be less pronounced. There is a large- 
scale decrease in flux after the instability and all dust disappears 
when the last planetesimal is ejected after ~ 100 Myr. 

The giant and terrestrial planet evolution differed only 
slightly between the mixed and seeds simulations. (Recall that 



the giant planets in the seeds simulations are identical to the 
mixed simulations in terms of their total mass and mass dis- 
tribution.) However, the giant planet instabilities in the seeds 
simulations appeared slightly weaker. The mean innermost giant 
planet eccentricity was 0.21 for the seeds simulations compared 
with 0.27 for the mixed simulationsQ The cause of this differ- 
ence appears to be late encounters between a giant planet and a 
seed embryo, after the end of giant planet-giant planet scattering. 
To test this, we sub-divide the 39 unstable bigseed simulations 
into the 22 simulations in which a giant planet underwent at least 
one planet-seed scattering event after the completion of planet- 
planet encounters and the 17 simulations that did not (i.e., the 
last planet-planet scattering event in these 17 simulations came 
after the last planet-seed scattering event). The simulations with 
late planet-seed encounters had systematically shorter instabili- 
ties (as measured by the time between the first and last planet- 
planet scattering events) and lower final giant planet eccentrici- 
ties. We think that what happens is that in these systems a planet- 
seed scattering event can lead to a small readjustment in one 
giant planet's orbit that separates it sufficiently from other giant 
planets to stabilize the system. However, we note that these sam- 
ples are relatively small and we cannot rule out that the weaker 
seeds instabilities are a product of small number statistics. 

Given the weaker giant planet instabilities in the seeds sim- 
ulations, terrestrial planet formation was correspondingly more 
efficient: only about 1/4 of the unstable seeds simulations de- 
stroyed all of their terrestrial material compared with more than 
40% for the unstable mixed simulations. However, the differ- 
ences between the planetary evolution in the seeds and mixed 
simulations are small compared with those between some of the 
other sets of simulations such as the lowmass and equal runs. 

Despite their influence on the outer planetesimal disk, seeds 
underwent little accretion. Among all fifty bigseed simulations, 
no seed accreted more than five planetesimals, and there was 
just a single seed-seed collision and 5 giant planet-seed colli- 
sions. There was slightly more accretion among the seeds in 
the smallseed simulations, which had a comparable rate of 
planetesimal-seed impacts but a higher rate of giant impacts, 
with 4 seed-seed collisions and 12 giant planet-seed collisions 
among the fifty simulations, although we note that all but 3 of the 
giant planet-seed collisions occurred very early and were proba- 
bly caused by the seed being placed on an orbit that was initially 
very close to a giant planet. 

Figure [8] shows the correlations between the dust-to-stellar 
flux ratio F/F srar (l 0pm) after 1 Gyr, and either the innermost gi- 
ant planet eccentricity or the total terrestrial planet mass. We plot 
results for the smallseed, bigseed, and mixed simulations. 
The most important difference between the seeds and mixed 
simulations is that the seeds simulations produce less dust at 
late times, especially warm dust that is observable at wave- 



2 A somewhat higher fraction of the bigseed simulations were un- 
stable compared with the mixed and smallseed simulations. This dif- 
ference is only moderately statistically significant and might be due in 
part to a small glitch in our initial conditions for the seeds simulations: 
the inner edge of the outer planetesimal disk was always 4 Hill radii 
exterior to the outermost giant planet but the icy embryos were always 
initially between 10-20 AU. Thus, in many cases the innermost one to 
two seeds were in immediate dynamical contact with the outermost gi- 
ant planet. This preferentially occurred when the outermost planet was 
very massive (and hence on a more distant initial orbit) and this glitch 
does not appear to have contaminated our results. In fact, the median 
instability time was later for the bigseed simulations than the mixed 
simulations, which is the opposite of what would be expected if the in- 
stabilities were systematically driven by seed-giant planet interactions. 
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Fig. 7. Evolution of a simulation that includes 5 icy embryos (shown in grey) in the outer planetesimal disk. Left: Orbital eccentricity 
vs. semimajor axis of each body in the simulati on. The size scales w ith the mass 1 / 3 and the color corresponds to the water content, 
using initial values taken to Solar System data dRavmond et alj|2004|) and re-calculated during impacts by mass balance. The giant 
planets are shown as the large black bodies and are not on the same size scale. The two surviving terrestrial planets have a = 0.71 
and 1.27 AU, m — 1.3 and 1.58 M ffl , and Myr-averaged e = 0.09 and 0.17, respectively. Both terrestrials have substantial water 
contents accreted from hydrated asteroidal material. The two surviving giant planets have a =4.1 and 24.4 AU, m =1.27 and 0.45 
My, and Myr-averaged e =0.32 and 0.33, respectively. Top right: The spectral energy distribution of the dust during five simulation 
snapshots. The dashed line represents the stellar photosphere. Bottom right: The ratio o f the dust-to-stellar fl ux F/F star at 25jum as 
a function of time. The rough Spitzer observational limit is shown with the dashed line (Trilling et al. 2008). An animation of this 
simulation is available at http://www.obs.u-bordeauxl.fr/e3arths/raymond/scatterSED^eedl3.mpg 



lengths shorter than about lOOyurn. Among just the unstable sub- 
set of simulations, 53 of 96 mixed simulations (55.2%^-^°) had 
dust fluxes that would be detectable with S pitzer after 1 Gyr o f 
evolution, i.e., with F/F mr (70pm) > 0.55 milling et al1l2008h . 
In comparison, 14 of 29 (48.3%^g* % ) unstable smallseed and 
16 of 39 (= 41%+8 '%) bigseed simulations were detectable at 
70/im after 1 Gyr. Given the statistical error bars, the decreased 
detection rate compared with the mixed simulations is only a l<x 
result for the smallseed simulations and 2<x for the bigseed 
simulations. 

The effect of the seeds is more apparent when consider- 
ing the stable simulations. The stable systems remain detectable 
at 98% or higher rate for each of the mixed, smallseed and 
bigseed sets of simulations. However, the actual dust bright- 
ness decreases dramatically for simulations with seeds. The me- 
dian F/F star (70pm) at 1 Gyr was 26.2 for the stable mixed sim- 
ulations, 5.3 for smallseed, and 1.7 for bigseed. Given the 



scatter, the difference between the mixed and seeds simula- 
tions is significant at the 3<x level and the difference between 
the bigseed and smallseed is significant at 5<x. 

At 25pm the differences are even more striking. In un- 
stable systems, 12/96 (12.5%+^*) mixed simulations, 2/29 
(6.9%+™*) smallseed simulations, and 0/39 (0 +4 4% ) bigseed 
simulations were above the Spitze r detection thresho ld of 
F/F slar (25pm) > 0.054 after 1 Gyr (iTrifling et al.ll2008l) . This 
constitutes a 1 -2<x difference. Among the stable systems, 51/56 
(9 1.1 mixed simulations were detectable at 25pm after 

1 Gyr but not a single stable bigseed or smallseed was de- 
tectable (0 +7 3% for the combined seeds simulations). 

With Spitzer's detection limits, debris disks at 70pm vastly 
outnumber those at 24pm. Around stars older than 300 Myr the 
fre quency of 24um dust excesses was estimated at 2.S% + ^tm 
bvlCarpenteret all d2009l) and at 1.9% + 1.2% bv lGaspar etaLl 
(120091) . However, these estimates are based on just a handful 
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Fig. 8. Correlations with the dust flux after 1 Gyr of dynamical and collisional evolution for the seeds simulations, as compared 
with the mixed simulations. The top panels show F/F star at 70jum vs. the eccentricity of the innermost surviving giant planet (left) 
and the total mass in surviving terrestrial planets (right). The bottom panels show the same comparisons but at 25/im. The mixed 
simulations are in black, the smallseed simulations in green, and the bigseed simulations in red. Filled circles represent stable 
simulations and open circles unstable simulations. 



of detections (from more than 100 targets). Amon g the 24/im 
detec tions are severa l systems such as 77 Cor vi (Wvat t et al.l 
12005b and HD 69830 (iBeichman et al.l2005HBrvden et al.l2 006) 
that appear to contain dust at ~ 1 AU; this dust has b een in- 
terpreted as either being transient (IWvatt et al.l l2007aj) or due 
to a very peculiar outcome of planet formation (IWvatt et al.l 
120101) . Thus, the frequency of systems with 24/im excess due 
to collisional equilibrium processes is significantly smaller than 
the quoted values. The frequen cy of dust excesses at 70/im is 
16.4%+§-8* dTrilling et al.l2008l) . This is at least 5-8 times higher 
than at 24yum, and the removal of systems with potentially tran- 
sient dust can only cause this ratio to increase. 

The mixed simulations produce an overabundance of 25/im 
dust excesses. After 1 Gyr of evolution, the frequency of de- 
tectable dust at 25/im was 91%^* and 12.5%+^* (1 - a er- 
ror bars) for stable and unstable systems, respectively. At 70/mi 
the frequency of detectable dust was 98.2%^-g* for stable and 



55.2%^j| a for unstable systems. Thus, the ratio of the fraction 
of systems that were detectable after 1 Gyr at 70//m to 25/itn is 
1.08 for the stable mixed systems and 4.4 for unstable mixed 
systems. The higher ratio for the unstable systems is a simple 
consequence of the instability preferentially clearing out the in- 
ner portion of outer planetesimal disks and leaving behind the 
colder part of the disk that does not emit much flux at 25/im. 
Nonetheless, no combination of these ratios for the mixed sim- 
ulations can match the observed ratio of 5-10 after 1 Gyr of evo- 
lutional 



3 We note that the unstable mixed systems with detectable 25/mi flux 
are all quite close to the detection limit (Fig. [8j, and the fraction of un- 
stable systems that is detectable at 25/im drops drastically to 2/96 = 
2. l%*QgS after 3 Gyr of evolution, and the ratio between the detectable 
frequency at 70/im to 25/mi increases to 25.5, well within the range 
allowed by observations. However, what is lacking in the mixed simu- 
lations is the ability to account for stable systems without 25/im flux, as 
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Fig. 9. Left: The fraction of systems that would be detectable with Spitzer (with F/F star (10fim) > 0.55 after 1 Gyr of collisional 
and dynamical evolution) as a function of the eccentricity of the innermost giant p lanet e g for the mixed and combined seeds 
simulations. The error bars are based on binomial statistics ("see lBur gasser et al. 2003). This essentially represents a horizontal slice 
through the top left panel of Fig [8] Right: The fraction of systems with 0.5 M e or more in surviving terrestrial planets as a function 
of F/F star (10fxm) > 0.55 (1 Gyr) for the mixed and seeds simulations. Systems with F/F star < 10~ 2 are included in the bin at 
F/F star » 10~ 2 . The Spitzer detection limit is shown as the dashed line. This represents a vertical slice through the bottom right 
panel of Fig. [8] 



The seeds simulations may explain the dearth of 25^m dust 
excesses because there is a natural suppression of the 25^m flux 
for both stable and unstable systems. From the 44 stable and un- 
stable smallseed simulations, only 2 were detectable at 25/mi. 
None of the bigseed simulations were detectable at 25//m. The 
reason for the lack of 25^m flux is that the viscous-like spread- 
ing out of the outer planetesimal disk acts to both deplete the 
inner part of the disk by inward scattering and to push the outer 
part of the disk to ever colder temperatures. The net effect is the 
near complete frustration of the 25//m flux. However, the detec- 
tion rates at 70/mi are higher than 50% for both the smallseed 
and bigseed simulations, including both the stable and unstable 
systems. Thus, the existence of ~ Earth-mass seeds in outer plan- 
etesimal disks may provide a natural explanation for the very low 
frequency of 25/iin excesses compared with 70/im excesses. 

As seen in Fig. [8] the anti-correlation between giant planet 
eccentricity e g and debris disks still holds for the seeds simu- 
lations. For eccentricities larger than 0.1, the dust fluxes show 
a rapid decrease for the sets of simulations with and without 
seeds because the dynamics of unstable giant planets dominates 
the survival of planetesimals. However, for smaller eccentric- 
ities there is a clear segregation: the bigseed systems have 
the smallest dust fluxes, the mixed systems have the largest, 
and the smallseed are in the middle. In this realm the stirring 
up of outer planetesimal disks dominates the dust flux, and as 
we've seen before the seeds simulations create a lower-mass 
and colder planetesimal disk than the mixed systems, leading to 
lower dust fluxes in proportion to the seeds' mass (not number). 
In fact, many of the unstable bigseed systems with modestly 
eccentric giant planets (e g ~ 0.1) have dust fluxes as high as 
the stable systems. However, despite their differences, Figure [9] 



the fraction of stable systems with detectable 25/um flux remains very 
high, at 83.9%^;™. 



shows that the fraction of systems that is detectable at 70^m as a 
function of e g is very similar for the seeds and mixed systems. 

The debris disk-terrestrial planet correlation still holds for 
the seeds simulations. As with the e g , the correlation between 
the dust flux at 70/im after 1 Gyr and the total surviving terres- 
trial planet mass is less evident because stable seeds systems 
that efficiently form terrestrial planets produce far less dust than 
their mixed counterparts. Again, for the stable systems there is 
a clear segregation between the sets of simulations, with the 
largest seed mass bigseed having the smallest flux. However, 
Fig.|9]shows that the fraction of systems that form at least 0.5 M e 
in terrestrial planets increases for both the seeds and mixed 
simulations. However, the seeds curve is 1 - 2<x higher than 
the mixed curve close to the detection limit. This reflects the fact 
that the seeds simulations deplete their outer planetesimal disks 
far more than the mixed simulations: only very calm systems 
preserve enough planetesimals to produce dust. The systems that 
are observed at a given dust flux therefore represent more stable 
systems for the seeds simulations than the mixed ones. Thus, 
the seeds simulations predict an even stronger correlation be- 
tween stars with observed debris disks and yet-to-be-discovered 
terrestrial planets. 

4.2. The width of the outer planetesimal disk (the widedisk 
simulations) 

The widedisk simulations allow us to test the effects of a 
higher-mass, wider outer planetesimal disk. Compared with the 
mixed simulations, the outer planetesimal disk in the widedisk 
simulations was twice as wide, 20 AU rather than 10 AU, and 
contained twice the total mass in planetesimals, 100M ffi instead 
of 50M ffi . The inner 10 AU of the outer planetesimal disk is 
therefore the same for the widedisk and mixed simulations (al- 
though the numerical resolution is halved in the widedisk runs) 
but the widedisk systems contain an additional 10 AU of plan- 
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Fig. 11. Left: The dust-to-stellar flux ratio at 1 Gyr at 70jum vs, the innermost giant planet's eccentricity for the widedisk (grey) and 
mixed (black) simulations. Filled circles represent stable simulations and open circles unstable ones. The Solar System is shown 
with the grey star. Right: Histogram of the fraction of systems that are detectable at 70fim as a function of the innermost giant 
planet's eccentricity; this is essentially a horizontal slice through the left panel. 
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Fig. 10. Orbital evolution of a widedisk simulation that under- 
went a Nice model-like instability. The evolution of the semi- 
major axes, perihelion and aphelion distances of the three gi- 
ant planets are shown in black, and the semimajor axes of the 
19 surviving outer planetesimals in grey (note that one addi- 
tional planetesimal survived with a = 59 AU). Two terrestrial 
planets formed in this system, at 0.71 AU (1.4M e )and 1.76 AU 
(0.13 Me): accretion was perturbed by the perturbations moder- 
ately eccentric inner giant planets at 4.7 AU (0.78My, e oscillates 
between 0.03 and 0. 17) and 7.8 AU (2AM j, e oscillates between 
0.06 and 0.1). 

etesimals. In addition, the outer boundary of each simulation - 
the limit beyond which particles are considered to be ejected - 
was 1 000 AU rather than 1 00 AU. 

The widedisk simulations behaved similarly to the mixed 
simulations in most respects. However, the more massive outer 
planetesimal disk did cause a few notable differences between 
the widedisk and mixed simulations. Given the much larger 



angular momentum reservoir contained in the outer planetesimal 
disk, ~Saturn-mass giant planets on eccentric orbits can be cap- 
tured in the outer system by planetesimal scattering. Figure [TUl 
shows the evolution of one such system, in which the two inner 
planets (M ilmer = 0.78Mj,M outer = 2AMj) startedjust interior to 
the 2:1 mean motion resonance. The giant planets' clearing out 
of the inner portion of the outer planetesimal disk drove the two 
inner giant planets across the 2:1 MMR after ~ 0.1 Myr. This 
caused a perturbative increase in their eccentricities and a cor- 
responding increase in the eccentricity of the outer, Saturn-mass 
planet. The outer planet's semimajor axis increased quickly and 
over the next 10 Myr its eccentricity slowly decreased by secular 
friction with the outer disk of planetesimals. Despite the pertur- 
bative evolution of the system, 19 planetesimals totaling 1.9 M e 
survived in the outer planetesimal disk (their orbital evolution is 
shown in grey in Fig. [TOt . Most of these started the simulation 
in the outer parts of the planetesimal disk (one at 18 AU). As 
the outer giant planet migrated outward, it shepherded many of 
these planetesimals in its 3:2, 2:1, and even its 3:1 mean motion 
resonances, located at 30.8, 37.2 and 48.9 AU at the end of the 
simulation, respectively. This is analogous to the shepherding of 
Kuiper belt objects such as Pluto durin g Neptune's planetesimal- 
driven migration (iLevison et al. 1 120081) . The surviving planetes- 
imal disk in the simulation from Fig.[l0]is massive enough to 
remain detectable at 70/im for 3 Gyr but is cold enough not to 
be detectable at shorter wavelengths. Planetesimal-driven migra- 
tion of a massive planet therefore represents another mechanism 

- in addition to the presence of seeds in outer planetesimal disks 

- to deplete outer planetesimal disks and to push them outward 
to colder temperatures. However, this mechanism operated effi- 
ciently in only a small fraction of widedisk simulations. 

Despite having identical giant planet initial conditions, a sig- 
nificantly smaller fraction of widedisk simulations went un- 
stable compared with the mixed simulations (40.7%^^ for 
widedisk vs. 63.2%*?-j? for mixed). This is because of the 
larger angular momentum reservoir in the outer planetesimal 
disks. For cases when the instability starts in the outer portion of 
the planetary system, as the outer planet's orbit becomes eccen- 
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Fig. 12. Left: The dust-to-stellar flux ratio at 1 Gyr at 10/j.m vs, the total terrestrial planet mass for the widedisk (grey) and 
mixed (black) simulations. Filled circles represent stable simulations and open circles unstable ones. The Solar System is shown 
with the grey star. Right: Histogram of the fraction of systems that contain 0.5 M ffl or more in terrestrial planets as a function of 
F I F star (70jjm) after 1 Gyr; this is essentially a vertical slice through the left panel. 



trie the disk's ability to transfer orbital angular momentum (as 
well as energy) to damp the outer planet's eccentricity and also 
to increase its semimajor axis increases for a more massive plan- 
etesimal disk. Thus, if all systems of giant planets formed on or- 
bits that would be unstable in the absence of planetesimal disks, 
their long-term orbital evolution should vary with the outer disk 
mass. The distribution of outer disk masses may therefore play a 
critical role in shaping the dynamics of inner planetary systems. 

FigureQ~T]compares the debris disk - giant planet eccentricity 
anti-correlation for the widedisk and mixed simulations. Both 
distributions follow the same shape as the other sets of simula- 
tions: at e g < 0.05 the vast majority of systems were dynamically 
stable and therefore have high fluxes, and for e g > 0.05 unstable 
systems dominate and the probability of having a significant dust 
flux decreases strongly with increasing e g . For stable systems, 
the 70/mi fluxes are generally 2-3 times larger for the widedisk 
systems because their outer planetesimal disks are more mas- 
sive and the dust flux is dominated by the outer regions of the 
disk where the collisional evolution is slow. However, the ratio 
of the median F / F star (70/j.m) after 1 Gyr for the widedisk to 
the mixed simulations was 2.7, but at 25yi/m the value was only 
1 .6. This shows that in the inner regions of the outer planetes- 
imal disk, where collisional evolution is faster, the higher-mass 
widedi sk disks approach the dust production level of the lower- 
mass mixed disks. 

The distribution of fluxes for the unstable systems is the 
same for the widedisk and mixed systems, meaning that the 
rate of survival of planetesimal disks is dominated by the giant 
planet dynamics rather than the initial conditions, even the total 
mass and radial extent. For these unstable systems, the distri- 
butions of instability times for the mixed and widedisk sim- 
ulations were indistinguishable. The distribution of the fraction 
of systems that is detectable with Spitzer is almost identical for 
the two sets of simulations: in both cases there is a plateau at 
~ 100% for e g < 0.05 and a sharp decrease toward higher ec- 
centricities. 

Figure[T2]compares the debris disk - terrestrial planet corre- 
lation for the widedisk and mixed simulations. Again, the two 



distributions have the same shape but the dust fluxes for the sta- 
ble systems (i.e., for those with the highest surviving terrestrial 
planet mass) are 2-3 times higher for the widedisk systems. The 
distributions are almost indistinguishable for systems with less 
than about 2 M ffi in surviving terrestrial planets. For both sets 
of simulations the fraction of systems with 0.5 M e or more in 
terrestrial planets increases monotonically with F/F star (lQ/j.m) 
at 1 Gyr. However, above the detection limit a slightly higher 
fraction of mixed systems contain terrestrial planets compared 
with widedisk. This is the opposite of the effect that we saw 
in section 4.1 for the seeds systems. The widedisk systems 
form very bright dust disks and require somewhat stronger gi- 
ant planet perturbations to decrease the flux below the detection 
limit. Thus, systems with bright debris disks at lOfim are less 
sensitive to the presence of terrestrial planets than the mixed 
systems. Note that this difference comes from the slightly higher 
fraction of widedisk systems with eccentric giant planets that 
produce 70/im excesses (see Fig. fTTb 

5. Effect of the gas disk during instabilities (the gas 
simulations) 

In the gas set of si mulations we adde d additional forces to the 
Mercury integrator ( Chamb ers! 1 1 9991) that acted on planetesi- 
mal and embryo particles to mimic the effects of the dissipat- 
ing gaseous protoplanetary disk from which the planets formed. 
In these simulations we included two effects: 1) aerodynamic 
gas drag due to the headwind felt by bodies orbiting at the 
Keplerian speed while gas orbits slower due to pressure sup- 
port. Aerodynamic drag acts most strongly on small objects, i.e. 
planetesimals, and leads to a rapid decay in eccentricity and in- 
clination as well as a slower decay of the semimajor axis; and 
2) tidal damping (also called "type 1 damping") due to gravita- 
tional interactions between objects and the disk. Type 1 damping 
increases for more massive bodies because it is caused by waves 
excited in the disk, meaning that this was an important source of 
dissipation for embryos but not for planetesi mals. Aerodynamic 
drag was calculated using standard models (Ad achi et al.lll976h 
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assuming planetesimals to be spheres with radii of 10 km. Type 1 
damping was included based on lin ear calculations for pl anets 
embedded within isothermal disks (Tanaka & Ward 2004). We 
included additional terms for large eccentricities and inclina- 
tions that were derived by ICresswell & Nel son (2008). We in- 
cluded type 1 damping but not type 1 migration both to maintain 
a clearer comparison with the simulations without gas forces and 
because eccentricity damping from the d isk is roughly 2 orders 
of magnitude faster than radial migration dTanaka & Wardl2 004; 
ICressweU & Nelsonll2008l) . 

We assumed the presence of an underlying gas disk 
that corres ponds to rough l y half the minimum -mass so- 
lar nebula dWeiden schilling Il977t lHavashil [l9 8 1 ), with sur- 
face density profile S(r) = 875 (r/lAlTT 3/2 gcmT 2 and 
vertical density distribution p(z) = exp(z 2 /zl), where 



ZnO) = 0.0472 (r/lAU) 5/4 AU (see also iThommes et al J 120031: 
iRavmond etall l2006at (Mandell et ail |2007l) . We note that the 
distribution of solid mass in our initial conditions is distributed 
according to a different density profile, 2 oc r _1 , in both the ter- 
restrial and outer planetesimal zones. We used a steeper radial 
surface density profile for the gas in order to increase the gas 
density in the inner disk to maximize the effect of gas drag on 
the survival of terrestrial bodies. However, we note that our ini- 
tial conditions for the inner and outer disks - which were chosen 
as represent approximate guesses for the Solar System's primor- 
dial disk - do not even follow the same global surface density 
profile because there is far too little mass in the terrestrial zone. 
The solution to this problem may lie with variations in the effi- 
ciency of planetesimal fo rmation at differe nt orbital radii within 
protoplanetary disks (e.g. lCha mbers 2010h. 

To model the final stage in the lifetime of the gaseous 
disk, the disk's surface density was dissipated linearly and 
uniformly in 500,000 years. This is slightl y longer than 
most estimates of the final dis sipation phase ( Simon & Prato 
19951 IWolk & Walter! [1995 IChiang & Murrav-Clavl l2007t 
Currie et all [2009) and should thus maximize the impor- 
tance of the gas disk phase. This situation is roughly con- 
sistent with models for dynamical in stabilities among plan- 
ets in the presence of gas disks dCha tteriee et al. 2008; 
Moeckel et alJ2008tlMatsumura et al.l2010l:lMarzari et al.l2 010: 
Moeckel & Armitagd 120121) . which predict that instabilities 



should preferentially occur late in the disk phase. The power- 
law gas density profile probab ly overestimates the amo unt of gas 
interior to the giant planets dCrida & Morbi delli 2007), so these 
simulations should provide an upper limit to the effects of gas on 
terrestrial bodies. The initial conditions for the gas simulations 
were drawn directly from the mixed set of simulations with one 
important change: the middle giant planet's eccentricity was in- 
creased to make its orbit cross the orbit of the innermost planet. 
This was to ensure that the system would be immediately unsta- 
ble so that we could test the effects of the relatively short-lived 
gaseous disk. 

The goal of the gas simulations is to test the effect of damp- 
ing from the gas disk on the dynamics and survival of rocky and 
icy bodies in the inner and outer planetary system. To accom- 
plish this, we want the giant planet instabilities to be the same as 
for the mixed set, to isolate the effects of the disk. Thus, we ne- 
glected type 1 and type 2 radial migration of giant planets (and 
embryos) in the gas simulations, although th ey would certainly 
occur in a re alistic minimum-mass disk (e.g jLin & Papaloizoul 
119861: IWardlfl997h . Indeed, in a more self-consistent setting the 
giant planets would likely be trapped in resonance at early times 
and the in stability would be t riggered by either ecce ntricity 
excitation dMarzari et all 120101: Libert & Tsiganisll201 ll) or the 
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Fig. 13. Total mass in surviving terrestrial planets as a function 
of the minimum perihelion distance of any giant planet during 
the simulation for the unstable mixed (black dots) and the gas 
(grey dots) simulations. 



dispersal of the gas disk ( Moeckel et al.ll200"8T: IChatteriee et al.l 
l2008HMoeckel & Armitagdl2012l) . 

The giant planets behaved similarly in the gas and the unsta- 
ble mixed simulations (note that in this section we compare with 
only the unstable portion of the mixed simulations because all 
of the gas simulations were unstable). The median eccentricity 
of the surviving giant planets in the gas simulations is slightly 
higher - 0.26 vs. 0.22. The reason for the stronger instabilities in 
the gas simulations is that they were initially placed on strongly 
unstable, crossing orbits and thus were unable to undergo weaker 
instabilities that can occur when they develop more slowly (and 
recall that no damping was felt by the giant planets in the gas 
simulations). The gas giants provided a comparable fit to the 
exoplanet distribution (p value from K-S test of 0.67 for gas, 
0.49 for mixed). A notable difference between the distributions 
is a population of low-eccentricity (e < 0.05) giant planets that 
is significantly more abundant in the unstable mixed than the 
gas simulations. This population was generated by weakly un- 
stable systems and these systems are the most efficient at both 
forming terrestrial planets and producing debris disks. Although 
the instabilities in the gas simulations occurred very early by 
design, the distributions of the duration of instabilities were vir- 
tually identical for the gas and unstable mixed simulations, ex- 
tending from 10 4 to 10 7 years with a median of slightly more 
than 300,000 years. This is due in part to the fact that we have 
not included appropriate drag forces acting on the giant planets 
as these are di fficult to estimate without hydrodynamical simu- 
lations (see e.g jMoeckel et al.ll2008T:lMarzari et al.ll2010h . 

The effect of the giant planets on terrestrial planet formation 
was similar for the gas and mixed simulations. Figure[T3lshows 
that the sculpting of the terrestrial zone by the giant planets - as 
measured by the minimum giant planet perihelion distance dur- 
ing the simulation - is the same for the unstable mixed and the 
gas simulations. The only slight difference comes from two gas 
simulations in which a giant planet came closer than 1 AU to 
the star (in one case for a prolonged period of almost 1 Myr) but 
that succeeded in forming a terrestrial planet. In both of these 
cases the surviving planet was roughly a Mars mass (although 
in one case the planet accreted another embryo) and underwent 
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Fig. 14. Distribution of the number of surviving terrestrial plan- 
ets per system for the unstable mixed (grey) and gas (dashed 
line) simulations. 



large-scale oscillations in eccentricity and inclination. The sur- 
vival of these planets may be due in part to gas drag from the 
disk, although they may simply represent a tail of the mixed dis- 
tribution. 

Figure [14] shows that the distributions of the number of sur- 
viving terrestrial planets in the gas and the unstable mixed sys- 
tems are very close (lcr is 5-8% for most bins as calculated us- 
ing binomial statistics; see Burga sser et al.l l2003). The gas and 
unstable mixed simulations each formed a mean of 1.2 planets 
planets per unstable system, although there is a slight tendency 
for more zero- and one-planet systems in the gas simulations. 
We attribute this to the fact that in the gas simulations the gi- 
ant planets started on orbits that were already strongly unstable 
and the instability involved the innermost giant planet with the 
strongest influence on the terrestrial planet zone. Thus, outward- 
directed instabilities and weakly unstable systems were less fre- 
quent in the gas simulations. 

Although the number of surviving planets was similar, the 
surviving terrestrial planets in the unstable mixed simulations 
were significantly more massive than in the gas simulations. 
The median terrestrial planet mass was 0.43 M ffi for gas and 
0.73 M ffi for mixed (counting only simulations that were inte- 
grated for >100 Myr and planets > 0.1 M e ). In addition, 28 of 
the 90 unstable mixed terrestrial planets were more than 1 M e 
(31.1%^*) compared with 4 of 33 (12.1%+™%) f or the gas 
terrestrial planets. The larger masses of the unstable mixed sim- 
ulations come from the contribution from weakly unstable sys- 
tems, i.e. those with minimum giant planet perihelion distances 
larger than about 4 AU in Fig. [13] For systems where a giant 
planet entered within 4 AU of the star, the two sets had the same 
median terrestrial planet mass. 

A small fraction of unstable systems produced asteroid belts 
without terrestrial planets. In these systems a number of rocky 
planetesimals were the only survivors in the inner planetary sys- 
tem, as all terrestrial embryos had been destroyed. This occurred 
in 14 of 299 (4.7 ± 1.2%) of unstable simulations across all the 
sets of simulations (excluding the lowmass simulations). The 
gas simulations had a slightly higher rate of production of as- 
teroid belt-only systems (3/45 = 6.7%+JqJ), presumably be- 
cause in a few cases gas drag was able to stabilize the orbits 



Fig. 15. The fraction of systems that contain > 0.5 M e in sur- 
viving terrestrial planets as a function of F/F slar (10fim) after 1 
Gyr for the unstable mixed (black) and gas (grey) sets of simu- 
lations. Error bars are lcr values calculated with binomial statis- 
tics. 



of planetesimals that were marginally unstable. These asteroid 
belts are low-mass, containing only 1-20 asteroid particles (each 
5 x 10~ 3 M e ), albeit typically on excited orbits with moderate 
eccentricities (e ~ 0.2 - 0.3) and inclinations (z ~ 20 - 30°). 
Given their rapid collisional evolution, these belts probably be- 
come quickly dominate by a few relatively large objects and are 
probably not detectable with current instruments. 

The surviving terrestrial planets in the gas and the unstable 
mixed simulations had similar median eccentricities e and incli- 
nations i (median values of e w 0.1 and i « 5 — 6°). However, 
the gas simulations tend to have significantly higher oscillation 
amplitudes in e and i. Although the median oscillation ampli- 
tudes are relatively close (median peak-to-peak e osc = 0.13 and 
i osc = 7° for gas vs. 0.10 and 5° for mixed), planets in the gas 
simulations are shifted to higher values. Again, this difference 
is simply due to the lack of weakly unstable systems in the gas 
simulations; when a cut in the minimum giant planet perihelion 
of <4 AU is applied the oscillation amplitudes are a match. 

The anti-correlation between the giant planet eccentricity 
and the dust flux at 70yum is very similar between the gas and un- 
stable mixed simulations. The positive correlation between the 
mass in surviving terrestrial planets and dust flux is also pre- 
served in the gas simulations. Figure [TBI shows the fraction of 
systems that formed at least 0.5 M 9 in terrestrial planets as a 
function of F/ F star (l Ofrni) at 1 Gyr. The two distributions are al- 
most identical. The only slight divergence is at small F/F star val- 
ues, where the mixed simulations are about lcr higher than the 
gas simulations. This is explained by the fact that the gas sim- 
ulations are inward-directed by design, because the instability is 
triggered by a close encounter between the inner two giant plan- 
ets (simply because our initial conditions put the middle giant 
planet on an orbit that crosses the inner one's). In contrast, the 
mixed instabilities include both inward- and outward-directed 
instabilities, i.e., instabilities that can be triggered in, and largely 
confined to, either the inner or outer parts of the system. Inward- 
directed instabilities that perturb the outer planetesimal disk 
are necessarily very strong, somewhat stronger than equivalent 
outward-directed instabilities that perturb the outer disk. This 
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only appears to be important (and then only at lcr) for the low- 
est values of F/F star , where the disk is most strongly depleted 
(note that the lowest bin includes systems with F/F star < 0.01). 
Thus, the presence of disk gas at the time of giant planet insta- 
bilities does not appear to have a significant effect on the debris 
disk-terrestrial planet correlation. 

We conclude that there are no strong systematic differences 
between the gas and the unstable mixed sets of simulations. 

6. Implications of dynamically active giant planets 

We now explore the implications of our results for expected cor- 
relations between extra-solar terrestrial planets, giant planets and 
debris disks. We emphasize that our conclusions are of course 
determined in part by our chosen initial conditions, which are 
poorly-constrained observationally. Several aspects of the initial 
conditions have considerable uncertainties: 1) the mass and mass 
distribution in the terrestrial planet zone, 2) the mass, mass dis- 
tribution and extent of the outer planetesimal disk, 3) the num- 
ber, masses and initial spacing of the giant planets, and 4) the 
relative masses and spacing of these three components. We have 
chosen what we consider to be reasonable values of the different 
components, but these values certainly vary from system to sys- 
tem and cover a far wider range than the subset included here. 
In addition, some systems with qualitatively different properties 
probably exist, such as disks with very widely-spaced giant plan- 
ets or in which planetesimals only form in narrow regions. We 
discuss the impact of these assumptions on our results in section 
6.4 below. 

With these limitations in mind, we now perform a simple ex- 
periment based on the results of our simulations. The goal of the 
experiment is to use the giant planets to match the observed exo- 
planet mass and eccentricity distributions and to then test differ- 
ent correlations within that framework. We construct two sam- 
ples of systems that provide an adequate match to observations 
using simple, non-fine-tuned mixes of different sets of simula- 
tions. We then explore the implications of those samples. 

6.1. Observational constraints 

Any sample we construct is constrained by observations of giant 
exoplanets, debris disks, and correlations between those two. We 
now summarize these key characteristics. 

First, we are constrained by the distribution of known gi- 
ant exoplanets, particularly those beyond 0.2 AU that have pre- 
sumably not been affected by tidal interactions with their host 
stars. The mass d istribution can be fit with a simple powe r law: 
dN/dM ~ M' L1 (lButler et alJl2006t lUdrv & Santosll2007l) . This 
is the mass distribution of the entire sample. It is likely similar 
to the ensemble-averaged mass distribution of planets prior to 
scattering (with some modification due to mass-dependent plan- 
etary ejections and collisions with the star) but, as discussed 
below, it need not be the mass distribution prior to scatter- 
ing in any individual system (i.e. on a system-by-system level, 
planet masses may be correlated). The frequency of giant plan- 
ets is very low (< 1%) close-in, but increases sharply at 0.5- 
1 A U and appears to remain at a high leve l out to at least 3 
AU dCumming et al.l2008L[Mavor et al.l201 ll) . The median giant 
exoplanet eccentricity is ~ . 25 and the distributio n extends to 
above 0.9 (Butl er et alJl2006HUdrv & Santosll2007l) . The eccen- 
tricity distribution is ind ependent of orbital distance (for plan- 
ets not affected by tides iFord & R asio 2008). In addition, ob- 
servations show that more massive giant planets (M p > Mj) 



have higher eccentricities than lower- mass giants (M D < Mj\ 
i Jones et alJl2005 iRibas & Miralda-Escud3l2007t IFord & Rasio 
l2008HWright et al.ll2009l) ~ 

Second, we are constrained by debris disk statistics, and 
we focus on observations at 70/im made primarily with Spitzer. 
Observations show that 16.4% + ^ of solar-type s tars have de- 
tectable dust emission at 70pm (Trillin g" et al.l2008h . There is no 
observed variation in this fraction with age, although the upper 
en velope of actual fluxes decreases for stars older than 1 Gyr or 
so dHillenbrand et alJl2008HCarpenter et"aT1l2009l) . 

Finally, we are constrained by any connection that might ex- 
ist between the presence of giant exoplanets and debris disks. 
Debris disks have been detected around more than 20 stars 
with known exoplanets. However, there is currently no cor- 
relation between the presence of planets and debris disks: 
the incidence of debris disks is about 1 5% for both stars 
with and without planets (see Table 1 in iKospal et all 12009: 
iMoro-Martfn et al.1 120071 IBrvden et all 120091) . We "also note 
that the strong observed correlation between the fraction of 
stars with currently-known exopl anets (hot Jupit e rs in partic- 
ular) and the stellar metaHici ty dGonzalezI 1 1997b ISantos et al.l 
1200 ll: iFischer & Valenti|l2005l) is not appa r ent in the sample of 
known debris disks (|Beichman et al. 2006t iGreaves et al. 2006; 
IBrvden et al.ll2006HK6spal et al.ll2009t) . 

6.2. Consistency with known giant planet properties 

The masses of individual giant planets, as well as the mass 
ratio between planets in a given system, are the most im- 
portant facto^g^venring_fhe outcome of planet-planet scat- 
tering (Raym ond et alj [2010). Equal-mass giant planets pro- 
vide the strongest instabilitie s for a given ma s s, and the most 
eccen tric surviving planets dFord et alJ 120031 : iRavmond et al.l 
2010). Scattered equal-mass planets are also more widely- 
spaced than p lanets with mass ratios of a few dRavmond et al.l 
2009a, 2010). And among equal-mass unstable systems, more 
massi ve planets yield large r eccentricities but smaller inclina- 
tions (Ray mond et alJl2010h . The dynamics of scattering is only 
weakly dependent on the planet masses; Neptune-mass planets at 
a few AU require far more close encounters to eject one another 
than Jupiter-mass planets but their final orbital distributions are 
similar. The number of giant planets also plays a role; in general, 
more giant pla nets lead to more scatter ing events and higher final 
eccentricities dJuric & Tre maine 20081). 

We construct two mixtures of our simulations to reproduce 
the observations: 

- Case A is based on the mixed simulations. The eccentric- 
ity distribution of surviving giant planets in the simulations 
(considering just the innermost planet as it provides the clos- 
est match to radial velocity observations) provides a quan- 
titative match to the observed distribution with a probabil- 
ity value p of 0.49 calculated from a K-S test. The best 
match is found by including only unstable systems, but the 
p value is still an acceptable 5-25% if the giant planet sam- 
ple includes a 5-10% contribution from stable systems, with 
a higher p for smaller contributions of stable systems (see 
Fig. 19 in Paper 1). Case A includes a 10% contribution from 
stable systems. Note that, since case A is built on the mixed 
simulations and that several other sets of simulations share 
the same giant planet characteristics (smallseed, bigseed, 
widedisk, and gas), variations on Case A can be con- 
structed by substituting a different set of simulations for the 
mixed set. 
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Fig. 16. The dust-to-stellar flux ratio at 70yum after 1 Gyr of evolution as a function of the eccentricity of the innermost giant planet 
e g j n for cases A and B. Systems in which the innermost giant planet is exterior to 8 AU have been excluded from this plot. 



- Case B is constructed from a combination of the equal and 
lowmass simulations. The simplest scenario to explain the 
observation that higher-mass planets have higher eccentrici- 
ties is for massive planets to form in system s with mul tiple, 
rough ly equal-mass planets (see section 5 of iRavmond et alj 
2010). At low planet masses the eccentricities are larger than 
observed, so to balance the sample a contribution of systems 
with lower-mass (M < Msm ~ Mj up ) planets with signifi- 
cant mass ratios is needed - these are represented with our 
lowmass set. The exact number of low-mass systems needed 
is poorly constrained. In practice, we divide the lowmass 
in two based on the mass of the innermost surviving giant 
planet (dividing at 50 M e ); case B includes an equal num- 
ber of systems at M < Mj up from the equal and unstable 
lowmass simulations. 

Cases A and B each provide a marginally acceptable match 
to the observed exoplanet eccentricity distribution. Case A 
matches the giant exoplanet distribution with p = 0.08 from K-S 
tests, and if we allow a 5-10 % increase in the num ber of planets 
with e — (as suggested by Zakamska et al.ll2010h Case B also 
matches the distribution, with p ^ 0.1. Gi ven the uncertainties in 
orbital fitting of exo planet eccentricities dShen & Turner| [2008; 
Zakamska et af1 l2010l) we do not attempt to fine-tune our sam- 
ples to better fit the observations. Case A naturally matches the 
observed mass distribution (except for a bias due to the mass de- 
pendence of ejected planets) and as the outcomes for the equal 
simulations were mostly mass-independent, weighting of differ- 
ent outcomes with the equal contribution is not necessary, and 
case B can be also consider ed to provide a mat ch to the mass 
distribution (see section 5 in Ray mond et al.ll20Tol for more de- 
tails). However, the two cases have different implications for 
the nature of planetary systems. In case A, all planetary sys- 
tems experience the same qualitative evolution because nearly 
all of them become dynamically unstable. In case B, the evo- 
lution of systems is divided according to the planetary masses: 
high-mass planetary systems undergo extremely violent insta- 
bilities, but the evolution of lower-mass systems is much calmer 
and many such systems are dynamically stable. Current observa- 



tions favor case B over case A becaus e it reproduces the higher 
eccentricities of more ma ssive planets (Ribas & Miralda-Escude 
120071; IWright et al.l l2009) while for case A high er-mass planets 
have lower eccentricities than low-mass planets (Raymond et al.l 
l2010h . 

6.3. Debris disk correlations in cases A and B 

The expected trends - an anti-correlation between giant planet 
eccentricity and debris disks and a positive correlation between 
terrestrial planets and debris disks - are clearly seen in both the 
case A and B systems. 

Figure [16] shows the dust flux at 70yum after 1 Gyr vs. the 
innermost giant planet's eccentricity for all the simulations in 
each case without the weighting described above (e.g., Case B 
still only contains the unstable lowmass simulations with in- 
ner planet masses greater than 50M ffl but there is no weight- 
ing between the equal and lowmass components). The anti- 
correlation is clear in the case A simulations but the large intrin- 
sic scatter in the lowmass component of case B makes the trend 
less evident for case B. The scatter is caused by the fact that the 
case B simulations are susceptible to planetesimal-driven radial 
migration that allows for a wide range in the depletion of the 
outer planetesimal disk depending on the outer planet's orbital 
history. 

The correlation between e g and F/F star (10pm) is clearer 
when the data are binned. Figure[T7](left panel) shows that cases 
A and B are very similar in terms of the fraction of systems that 
is detectable at 70pm as a function of e g : both cases show the ex- 
pected clear anti-correlation between debris disks and eccentric 
giant planets including a rapid decrease in the detectable fraction 
for e g > 0.03 - 0.1. When considering the fraction of systems 
with e g > 0.1 as a function of F / F star (70pm) (right panel of 
Fig.fTTl there is an offset of roughly 1 — cr between the two cases 
that is again due to the larger inherent scatter in the case B sim- 
ulations, in particular the high-mass, unstable component of the 
lowmass simulations. 
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Fig.17.Left: The fraction of systems that would be detectable with Spitzer (i.e., with F/ F star (70/jm) > 0.55 after 1 Gyr of collisional 
and dynamical evolution) as a function of the eccentricity of the inn ermost giant planet e g , for cases A (black) and B (grey). The 
error bars are based on binomial statistics (see[Burgasser et al. 2003). This essentially represents a horizontal slice through Fig[T6l 
Note that cases A and B include all relevant simulations; there is no weighting of stable/unstable or equal/lowmass simulations. 
Right: The fraction of systems with e g > 0.1 as a function of F / F slar {10nm) > 0.55 (1 Gyr), again for cases A and B. Systems with 
F IF star < 10~ 2 are included in the bin at F/F star ~ 10~ 2 . The Spitzer detection limit is shown as the dashed line. This represents a 
vertical slice through Fig. [16] 
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Fig. 18. The dust-to-stellar flux ratio at 70/itn after 1 Gyr of evolution as a function of the total mass in surviving terrestrial planets 
for cases A and B. Again, the Solar System is represented by the grey star. 



We now turn to the debris disk - terrestrial planet cor- 
relation shown in Paper 1. This correlation arises naturally 
from the destructive role of dynamically active giant planets 
for both terrestrial planet formation and the survival of outer 
planetesimal disks. Indeed, the giant planets' role in terrestrial 
planet formation is almost purely antagonistic: giant planets 
may quench or stifle terrestrial planet formation but have not 
been sho wn to help it along in any significant way (see Paper 1 
and also IChambers & Cassenl I2002T i Levison & Agnorl 1 2003 ; 
Ray mond et al.l l2004t iRavmondl l2006t iRavmond et al.l 12006a, 



2009c), except in some circumstanc es to pro mote runaway 
growth of planetesimals (Kortenkamp et al. 2001). Nonetheless, 
if we assume that giant planets generally form outside the snow 
line, moderately high eccentricities are needed before the impact 
on terrestrial planet formation becomes deleterious. 

Figure [18] shows the terrestrial planet-debris disk correla- 
tion for cases A and B. As before, the correlation between 
F I F star (70fim) after 1 Gyr and the total mass in surviving ter- 
restrial planets is clearer for case A. Again, this comes from the 
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low mass simulations' much larger scatter in F/F star and some- 
what lower typical values for F/F star . 

When the data are binned, the terrestrial planet-debris disk 
correlation is clearer for both cases. Figure[19](left panel) shows 
the fraction of systems with F '/ F SMr (70yum) above the Spitzer de- 
tection limit as a function of the total terrestrial planet mass (i.e., 
a horizontal slice through Fig. fTST ). The correlations are similar 
for cases A and B, with only small differences for small terres- 
trial planet masses, for which a higher fraction of systems re- 
mains detectable for case A. We attribute this difference to the 
more dynamic evolution of case B systems, which deplete or de- 
stroy their outer planetesimal disks at a much higher rate than 
case A systems. Cases A and B follow nearly identical trends 
in terms of the fraction of systems with at least 0.5 M e in ter- 
restrial planets as a function of F/F star (l Qfim) (i.e., a vertical 
slice through Fig.lT8b. The only modest discrepancy that at very 
low fluxes (F I F star (10im\) < 0.05) there are more case A sys- 
tems with terrestrial planets. Once again, we attribute this to the 
more dynamic evolution of case B systems: the typical case B 
system that destroys its outer planetesimal disk is more likely to 
also destroy its inner rocky material than a corresponding case 
A system because in this regime case B is dominated by the very 
violently unstable equal systems. 

To summarize, we conclude from Figs [T6TT81 that the debris 
disk - eccentric giant planet anti-correlation and the debris disk 
- terrestrial planet correlation are clear in both cases A and B. 

6.4. Discussion 

Could our two predicted correlations be artifacts of our initial 
conditions? Is there any reasonable scenario that could remove 
these correlations? 

The debris disk - eccentric giant planet anti-correlation ex- 
ists because giant planet instabilities tend to clear out outer plan- 
etesimal disks, mainly by dynamical ejection. For eccentric gi- 
ant planets not to be anti-correlated with debris disks, some- 
thing fundamental about our proposed scenario must change. 
To start with, other mechanisms have been proposed to explain 
the large eccentricities of the observed exoplanets (for an ex- 
haustive list, see Section 1 of Ford &~Rasiol 1200 8). However, 
to date the planet-planet scattering model is the only mecha- 
nism that has been shown to be physically viable and to fully 
reproduce the curre ntly-observed charac t eristics of the giant ex- 
oplanet population (Ford & Rasio 2008; Chatterjee et al. 2008; 
I Juric & Tremainell2008b IRavmond et al.H20 1 Oh . 

The initial conditions of our simulations certainly play a role 
in producing the correlations we find; for example, if there exist 
wide radial gaps between the giant planets and outer planetesi- 
mal disks then the giant planets' influence on debris disks would 
be weaker. The giant planets would appear to be basically irrel- 
evant for the existence of debris disks. However, even a distant 
giant planet has a destructive influence on a planetesimal disk by 
increasing eccentrici ties sufficiently that int er-particle collisions 
become destructive dMustill & Wv att 2009). A simple fit to the 
observed distribution of debris disks using a self-stirred model 
like the one presented in section 2 shows th at outer planetesi- 
mal d isks are typically located at 15-120 AU dKennedv & Wvatil 
1201(H) . Thus, it may be possible for debris disks and giant planets 
at a few AU or less to appear unco r rected (as is the ca se in the 
current sample iBrvden et alj 2009: K ospal et al .1120091) but the 
anti-correlation between eccentric giants and debris disks should 
be clear for giant planets beyond roughly 5-10 AU. In addition, 
there should be a minimal orbital radius (and a maximum tem- 
perature) for an outer planetesimal disk in a given system, de- 



pending on the giant planet architecture. The statistics of debris 
disks in known giant planet systems are currently too poor to 
provide constraints (see discussion in section 5.2 of Paper 1). 

The debris disk - terrestrial planet correlation exists for three 
reasons. First, terrestrial planet formation is less efficient in the 
presence of eccentric giant planets. Second, the destruction rate 
of outer planetesimal disks by ejection increases under the influ- 
ence of eccentric giant planets. Third, we have assumed that any 
system that can form terrestrial planets can also produce a debris 
disk. 

The fir st reason is well-established from many dynami- 
cal studies (IChambers & Cassenll2002l: iLevison & Agnor|[2 003: 
I Ravmond et alj 120041; iRavmondl 120061 IRavmond et all 12006a, 
2009c; Morish imaet alJl2010h . As discussed above, the dynam- 
ics behind the second reason is sound and the only reasonable 
way to negate or temper the eccentric giant planet - debris disk 
anti-correlation is for there to exist a large radial gap between gi- 
ant planets and outer planetesimal disks. The third reason is hard 
to constrain. As a blanket statement this assumption is almost 
certainly i ncorrect, as the observed frequency of d ebris disks 
of -16% dTrilhng et alJl2008t ICarpenter et alJl2009h is smaller 
than the cu rrently-estimated frequency of close-in super Earths 
of 20-50% dHoward et alJ20ldl201ltlMavor et alJl201 ll) . Thus, 
there are probably many systems that can form terrestrial plan- 
ets but without a massive outer planetesimal disk. The reason 
for the lower frequency of outer disks is unclear; it could be 
related to external pertur bations from passing s tars d uring the 
embedded cluster phase (Malmbe rg et alj 120071 l20lH) . Or per- 
haps outer planetesimal disks simply are not a common ini- 
tial condition for planet formation, although that assertion ap- 
pears to be at odds with the high freq uency of disks around 
young stars (e.g Hillenb rand et al.l[2008h . One can imagine that 
there could also exist systems with outer planetesimal disks 
but very little mass in the inner disk, for example if a giant 
planet migrated through and depleted the i nner disk (although 
Raymond et alJ l2006b: Man dell et alJl2007l shows that this de- 
pletion is much weaker than one would naively expect). Indeed, 
in our own Solar System Jupiter's inward-then-outward migra- 
tion may have remove d most of the mass from the asteroid 
belt dWalshet alJl20TTh . However, there is no evidence in the 
exoplanet distribution for systematic depletion of inner disks by 
giant planet migration. Indeed, the radial distributio n of giant 
exopl anets increases sharply beyond about 1 AU (Mayo r et alj 
201 1) such that terrestrial planets closer than 0.5-0.7 AU to their 
stars are probably only weakly affected by giant planets, at least 
those with low to moderate eccentricities. In contrast, inner disks 
certainly have a wide mass range and high-mass disks may form 
super Earths or Neptune-like planets rather th an Earth-like plan- 
ets dlkoma et alJl200U IRavmond et al.ll2008bl) . 

We think that the most likely interpretation of current obser- 
vations and theory is roughly as follows. Protoplanetary disks 
start with a variety of masses and mass distributions, and are 
subsequently divided into different regions by any giant plan- 
ets that form. Formation models suggest th at giant planets may 
preferentially form a t a few to ten AU d Kokubo & Ida 2002; 
Thommes et al. 200^ ILevison et al.ll2010l) . essentially dividing 
their disks into the distinct regions we have assumed: the inner 
terrestrial zone, the giant planet zone and the outer planetesimal 
disk. In some systems these zones are not cleanly separated, for 
example if the giant planets migrate very far outward or inward 
or if a relatively close stellar encounter disrupts the outer plan- 
etesimal disk. Indeed, observations suggest that there is probably 
an abundance of systems with inner terrestrial planet-forming 
disks but without the outer planetesimals to produce debris disks. 
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Fig. 19. Left: The fraction of systems that would be detectable with Sptizer (F/F star (10fjm) > 0.55 after 1 Gyr) as a function of the 
total mass in surviving terrestrial planets for cases A (black) and B (grey). This represents a horizontal slice through Fig[18] Right: 
The fraction of systems with 0.5 M ffl or more in terrestrial planets as a function of F/F star (70fj,wi) > 0.55 (1 Gyr) for cases A and B. 
Systems with F/F star < 10~ 2 are included in the bin at F/F slar w 10~ 2 . The Spitzer detection limit is shown as the dashed line. This 
represents a vertical slice through Fig. [19] 



However, we see no evidence or clear theory to contradict our 
assumption that all or at least most stars with debris disks also 
have inner disks of protoplanets. Thus, we think that debris disks 
can indeed act as signposts for systems that should have formed 
terrestrial planets. In contrast, in systems without debris disks 
eccentric giant planets may act as signposts of terrestrial planet 
destruction. 

To summarize, we see no compelling argument against the 
assumptions we have made and we think that our two proposed 
correlations are reasonable. 

7. Conclusions 

In Paper 1 dRavmond et alj|201 ll) we showed that old solar-type 
stars with bright cold dust correlate strongly with dynamically 
calm environments that are conducive to efficient terrestrial ac- 
cretion. The fact that both the inner and outer parts of plane- 
tary systems are sculpted by what lies in between - the giant 
planets - yields a natural connection between terrestrial planet 
and debris disks. We predicted two observational correlations: 
an anti-correlation between eccentric giant planets and debris 
disks, and a positive correlation between terrestrial planets and 
debris disks. We also showed that the Solar System appears to 
be a somewhat unusual case in terms of having a rich system 
of terrestrial planets but a severely depleted Kuiper belt with lit- 
tle cold dust, which is probably a result of the outward-directed, 
relatively wea k dynamical instability that caused the late heavy 
bombardment dMorbidelli et aLll2010h . 

In this paper we used six new sets of simulations to test the 
effect of several system parameters on the results from Paper 1 . 
None of the parameters qualitatively changed the eccentric giant 
planet-debris disk anti-correlation or the debris disk-terrestrial 
planet correlation but their effects were diverse and interesting: 

- Low-mass giant planets undergo planetesimal-driven migra- 
tion that radially spreads the giant planets (the lowmass sim- 
ulations). Systems with low-mass giant planet are very ef- 
ficient at forming terrestrial planets and also at producing 



abundant observable dust simply because their weak giant 
planet perturbations do not stifle these processes. Systems 
with low-mass outer planets often produce isolated belts of 
planetesimals orbiting between the giant planets. The prob- 
ability of a system containing such a belt increases strongly 
with decreasing outer giant planet mass and increasing sepa- 
ration between the giant planets. 

- In contrast, systems with equal-mass giant planets undergo 
the strongest instabilities. In more than 2/3 of the equal 
systems, all terrestrial material was destroyed, and all plan- 
etesimal disks were ejected in a similar fraction of cases 
(though not with a one to one correspondence). Given the 
tendency f or more massive ex oplanets to have more eccen- 
tric orbits dWright et al J 12009). the equal systems may be 
representative of the high-mass (M p > Mj) exoplanet sys- 
tems. 

- The presence of 0.5 - 2M e objects in outer planetesimal 
disks cause the disks to radially spread (the seeds simu- 
lations). The inward-spreading part of the disk is removed 
by interactions with giant planets (or, in systems without gi- 
ant planets, presumably by accelerated collisional evolution) 
whereas the outer part of the disk spreads to larger orbital 
distances. Such cold disks are detectable at long but not short 
wavelengths. The presence of seeds may therefore explain 
the very low frequency of stars with 25/ um excess compared 
with the frequency of 70/^m excesses dBrvden et all [2006: 
iTrilling et al 12008^ ICaipenter et al.ll2009l) . 

- The widedisk simulations showed that a more massive and 
extended outer planetesimal disk stabilizes a significant frac- 
tion of systems because planet-planetesimal interactions al- 
low higher-mass unstable giant planet systems to radially 
spread out and damp their eccentricities. Thus, if all systems 
of giant planet systems form in similar, near-unstable con- 
figurations, the outer planetesimal disk mass may be a key 
factor in the fraction of systems that end up being unstable. 
More massive outer disks produce larger quantities of dust, 
in particular at long wavelengths that are dominated by the 
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outermost part of the planetesimal disks where the collisional 
timescales are the longest. 
- The presence of disk gas for the first 0.5 Myr of the sys- 
tem's evolution does not have a strong effect on the outcome 
(the gas simulations). The debris disk-terrestrial planet cor- 
relation and the eccentric giant planet - debris disk anti- 
correlation are not affected. 

We constructed two samples of simulations called cases A 
and B that matched the observed exoplanet mass and eccentric- 
ity distributions. Case A was built on the mixed simulations and 
case B on a combination of the equal and the high-mass un- 
stable component of the lowmass simulations. These cases at- 
tempt to bracket the likely initial conditi ons for planet- planet 
scattering in the known systems (e.g., IChatteriee et"aT 



et-planet 

all 120081: 



Uuric & Tremaiiiel [2008: Ra ymond et all 1201 Ol) . Each of these 
cases clearly show the eccentric giant planet - debris disk anti- 
correlation and the debris disk-terrestrial planet correlation. 
The only plausible alternative to the eccentric giant planet 

- debris disk anti-correlation is for outer planetesimal disks to 
be located far away from the giant planets and thus to be ex- 
tremely cold (see section 6). Given that the known exoplanets 
are dominated by relatively close-in planets (typically within a 
few AU) and that the planete simals responsible for t he known 
debris disks are at 15-120 AU dKennedv & Wvattll2010h . it is not 
surprising that the expected eccentric giant planet - debris disk 
anti-correlation has not yet been detected (iBrvden et al-l feOOS*. 
see also section 5.2 of Paper 1). However, we predict that it will 
be seen in upcoming datasets that include more distant planets. 

The only plausible alternative to the debris disk-terrestrial 
planet correlation is for star with debris disks to be systemat- 
ically depleted in their inner (terrestrial) regions and for outer 
planetesimal disks to be systematically cold enough to be uncor- 
rected with eccentric giant planets. Although it is easy to imag- 
ine scenarios in which this may occur - for example, a system in 
which an inward-migrating giant planet depletes the inner disk 

- there is no observational evidence or self-consistent theoreti- 
cal model that suggests that such scenarios should be common 
or correlated with debris disks. Thus, we predict that upcoming 
datasets will find a strong correlation between debris disks and 
the presence of terrestrial planets. 

The frequency of de bris disks around S un-like stars older 
than 1 Gyr is -16% dTrilling et all l2008h . which is lower 
tha n the observed frequency of close-in super Earth plan- 
ets dHoward et alJl2010L 1201 lb iMavor et al.ll201 ll) . Thus, there 
is certainly a large population of systems that can form terres- 
trial planets but without the outer planetesimal disks needed to 
produce debris disks. In these systems, direct radial velocity ob- 
servations or constraints on the giant planet architecture may of- 
fer the best insight for constraining the existence of terrestrial 
planets. But stars with debris disks most likely also had abun- 
dant material for building terrestrial planets. Thus, debris disks 
can indeed act as signposts of (past) terrestrial planet formation. 

To sum up, our main result is a prediction that debris disks 
should be anti-correlated with systems containing eccentric gi- 
ant planets and correlated with the presence of terrestrial planets. 
Solar-type stars with bright cold debris disks and no giant plan- 
ets are excellent candidates to search for Earth-like planets. In 
contrast, systems without debris disks and with eccentric giant 
planets are probably not good candidates for terrestrial planets. 
Upcoming observations will test our predictions. 
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